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1 .  INTRODUCTION 


This  report  describes  the  modeling  activities  carried  out  by  Physical 
Sciences  Inc.  (PSI)  during  the  first  year  of  our  ONR  Contract  No.  0014-87-C- 
0161.  The  objective  of  this  contract  is  to  develop  models  which  can  ulti¬ 
mately  be  used  to  predict  the  thermal  response,  ablation,  and  fragmentation  of 
"hard"  and  "soft"  biological  tissue  exposed  to  free  electron  laser  (FEL) 
radiation.  Most  of  our  efforts  thus  far  have  been  focused  on  the  "hard" 
tissue  problem,  i.e.,  kidney  stones,  gall  stones,  calcified  plaque.  A  summary 
of  the  principal  ablation  mechanisms  under  consideration  is  presented  in 
Figure  1-1. 

The  models  developed  are  being  used  to  identify  the  dominant  physical 
processes  and  key  parameters  controlling  the  ablation/fragmentation  process. 
Further,  to  validate  the  models,  we  are  comparing  their  predictions,  whenever 
possible,  to  available  laser/tissue  interaction  data  that  have  been  generated 
in  the  medical  research  community.  Model/data  comparisons  are  vital  for 
guiding  new  model  developments.  Similarly,  the  models  help  point  out  the  need 
for  additional  critical  experiments. 
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During  this  first  year  of  the  program,  modeling  efforts  have  been 
initiated  in  four  principal  areas:  (1)  further  development  and  application  of 
1-D  and  2-D  numerical  heat  conduction  codes  to  model  the  thermal  damage  and 
ablation  of  soft  tissue  exposed  to  infrared  laser  energy  (Section  2); 

(2)  modeling  of  plasma  initiation  mechanisms  in,  and  resulting  plasma  spectral 
emissions  from,  pulsed  laser-irradiated  kidney  stones  and  gallstones 
(Section  3);  (3)  modeling  of  .  aser-driven  shock-induced  fragmentation  of 
kidney  stones  and  gallstones  (Section  4);  and  (4)  modeling  of  the  above¬ 
surface  hydrodynamics,  i.e.,  bubble  growth  and  acoustic  wave  propagation 
resulting  from  plasma-mediated  laser  energy  deposition  above  a  hard  tissue 
surface  in  a  fluid  medium. 

Wherever  possible,  the  predictions  of  the  models  have  been  compared  to 
available  data.  The  only  exception  is  in  the  case  of  the  mathematical  models 
describing  the  above-surface  fluid  mechanics  during  laser  lithotripsy. 
Numerical  solutions  to  the  relevant  equations  were  not  available  in  time  to  be 
incorporated  into  this  report.  Preliminary  numerical  results  will,  however, 
be  presented  at  the  upcoming  contractors'  review  meeting. 
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2.  NUMERICAL  MODELING  OF  THE  LASER-INDUCED  THERMAL  RESPONSE 


AND  ABLATION  OF  BIOLOGICAL  TISSUE 

Ve  present  in  this  section  detailed  calculations  of  the  thermal  response 
and  ablation  of  biological  tissue  exposed  to  lOy  laser  radiation.  For  most  of 
the  calculations  shown,  tissue  response  was  modeled  using  generalized  computer 
codes  that  were  originally  developed  to  model  the  laser-induced  response  of  a 
variety  of  non-biological  materials.  The  models  are  applicable  at  any 
wavelength  for  which  the  absorption  cross  section  in  the  tissue  is  large 
compared  to  the  scattering  cross  section. 

In  this  initial  study,  we  have  modeled  soft  tissue  as  a  single  component 
material  with  thermal  and  optical  properties  very  similar  to  that  of  water. 

The  predominant  phenomena  treated  are  in-depth  laser  absorption,  thermal 
conduction,  and  surface  vaporization.  In-depth  thermal  damage  of  tissue  is 
calculated  by  incorporating  into  the  codes  a  reaction  for  the  thermal 
denaturation  of  tissue  protein. 

Model  calculations  are  presented  for  exposures  ranging  from  low  pover  C V 
exposures  of  a  few  watts  per  cm^  and  exposure  times  of  tens  of  seconds  to 
single-pulse  exposures  with  peak  irradiances  of  several  megawatts  per  cm^  and 
microsecond  exposure  durations.  Where  possible,  model  predictions  are  com¬ 
pared  to  available  experimental  data  and  the  implications  of  the  comparisons 
are  discussed.  Recommendations  for  future  modeling  improvements  are  also 
presented. 

2.1  Introduction 


In  recent  years,  there  have  been  a  significant  number  of  research  studies 
directed  at  understanding  the  response  of  biological  tissue  to  laser 
irradiation. I  A  major  area  of  interest  has  been  the  laser-induced  thermal 
response  and  ablation  of  tissue  for  surgical  treatments,^-^  with  perhaps  the 
most  prominent  applications  being  tissue  cutting  and  tissue  welding  proce¬ 
dures.  The  motivation  has  been  to  help  define  the  laser  treatment  control 
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parameters  for  best  achieving  the  desired  beneficial  effects  of  the  laser, 
i.e.,  controlled  heating,  coagulation,  cutting  and  hemostasis,  and  minimizing 
the  undesired  effects,  i.e.,  unwanted  damage  to  adjacent  healthy  tissue. 

We  report  here  the  results  of  efforts  ve  have  initiated  to  develop 
detailed  models  of  laser  interactions  with  biological  tissue.  The  overall 
objective  has  been  to  establish  a  well-founded  theoretical  understanding  and 
validated  predictive  capability  for  characterizing  a  wide  variety  of  medical 
treatments  involving  laser/tissue  interactions.  The  ultimate  aim  is  to 
develop  predictive  models  that  are  both  reliable  and  versatile  and  which  can 
be  used  for  optimizing  the  design  of  medical  lasers  as  well  as  the 
surgical/clinical  techniques  by  which  they  are  applied. 

The  specific  approach  has  been  to  apply  detailed  heat  conduction  and 
ablation  codes  originally  developed  for  use  with  non-biological  materials^  to 
the  problem  of  biological  tissue  interactions.  For  the  present  studies,  the 
scope  has  been  focused  on  modeling  primary  tissue  response  to  10ym  laser 
irradiation.  However,  the  models  developed  can  be  applied  at  any  laser 
wavelength  for  which  absorption  dominates  over  scattering  and  the  absorption 
depth  in  tissue  is  known.  For  example,  work  is  now  underway  at  modeling 
tissue  interactions  at  the  Erbium: YAG  laser  wavelength  (A  -  3  ym). 

Interactions  for  a  wide  range  of  laser  temporal  waveforms  and  intensities 
have  been  studied.  The  models  have  been  applied  to  single-pulse, 
repetitively-pulsed,  and  continuous  wave  (CW)  irradiations.  We  present  here 
only  the  results  for  CW  and  single-pulse  exposures.  Exposure  times  have  been 
investigated  from  the  single-pulse  microsecond  regime  to  CW  exposures  of  tens 
of  seconds.  Irradiances  have  been  varied  from  as  low  as  2  W/cm^  to  peak 
irradiances  as  high  as  3  MW/cm^.  Although  most  of  the  calculations  have  been 
performed  in  the  limit  of  one-dimensional  heat  conduction,  predictions  are 
also  presented  in  which  two-dimensional  heat  conduction  effects  have  been 
modeled  in  detail.  The  latter  were  obtained  using  a  three-dimensional  heat 
conduction  and  ablation  model. 

Tissue  material  properties  have  been  modeled  rather  simply.  Tissue  has 
been  treated  as  a  single-component  homogeneous  material  with  thermal  proper- 


2-2 


ties  quite  similar  to  those  of  vater.  This  assumption  is  very  reasonable  for 
high  water  content  tissue  such  as  fresh  dermal  or  mucosal  tissue. 

Vhere  possible,  model  predictions  have  been  compared  to  available  experi¬ 
mental  data.  The  most  detailed  comparisons  have  been  made  with  the 
experiments  of  Takata  et  al.^  and  Walsh  et  al.b  The  former  studies 
investigated  the  thermal  damage  of  pig  skin  induced  by  CV  CO2  laser  radiation 
at  low  power  irradiations  of  1  <  I  <  9  W/cm^.  In  the  latter  investigations, 
tissue  ablation  was  studied  using  a  pulsed  TEA  CO2  laser. 

2.2  Description  of  Models  Used 

With  the  exception  of  a  simple  analytical  laser  blow-off  model  (presented 
in  a  following  subsection),  all  the  modeling  performed  in  this  study  was 
carried  out  using  numerical  heat  conduction  and  ablation  codes**  developed 
previously  at  Physical  Sciences  Inc.  (PSI).  These  codes  have  been  exercised 
extensively  over  many  years  to  analyze  the  response  of  a  wide  range  of  non- 
biological  materials  exposed  to  a  variety  of  high  power  lasers.  Two  different 
codes  have  been  used:  1)  a  generalized  ablation  code  that  treats  heat 
conduction  in  only  one  dimension,  but  which  incorporates  the  ability  to  model 
in  detail  a  wide  range  of  ablation  phenomena;  and  2)  a  three-dimensional 
ablation  code  that  can  model  radial  as  well  as  axial  heat  losses,  but  is  more 
limited  in  its  treatment  of  material  ablation.  In  general,  the  single-pulse 
calculations  have  been  performed  using  only  the  one-dimensional  ablation  code. 
For  the  single-pulse  calculations,  a  one-dimensional  analysis  is  clearly 
appropriate  as  radial  heat  conduction  can  be  shown  to  be  negligible  over  the 
time  scale  of  an  individual  pulse. 

A  brief  description  of  the  numerical  codes  is  presented  below. 

2.2.1  One-Dimensional  Laser  Effects  Code 


This  code  is  a  transient,  thermal  conduction  model  that  can  be  used  to 
treat  a  system  comprised  of  up  to  three  components.  For  organic  composite 
materials,  for  example,  it  treats  fiber  (or  a  non-decomposing  homogeneous 
component),  resin  and  char.  The  latter  two  are  associated  with  the  organic 
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(pyrolyzing)  components  of  the  system,  but  have  not  been  employed  to  date  for 
medical  applications  since  ve  have  used  a  simple,  single-component  description 
of  tissue  for  our  initial  studies.  The  one-dimensional  code  is  complemented 
by  a  variety  of  subroutines  dealing  with  specific  laser-induced  phenomena  that 
one  may  wish  to  model  according  to  the  application. 

The  external  laser  heating  can  be  CV,  single  pulse,  or  repetitively 
pulsed.  Energy  transfer  mechanisms  that  can  be  treated  in  the  code  are: 

o  Surface  and/or  internal  absorption  of  laser  energy 
o  Surface  vaporization 

o  Internal  pyrolysis 

o  Surface  thermal  radiation 

o  Energy  transported  by  heated  gaseous  pyrolysis  products,  and 

o  Convective  heating  (cooling)  at  the  surface  in  external  cross  flow. 

The  conservation  and  transport  of  energy  is  governed  by  the  following  equation 

3T  3  f„3T)  .  n  .  „  .  3T  dl 

*31  *  35  lK  '  prQp  “  PfHV  +  mgCg^  '  3x  (1 

where  p  is  the  density,  c  is  the  specific  heat,  t  is  the  time  coordinate,  x  is 
the  space  coordinate;  K  is  the  thermal  conductivity,  pr  is  the  time  rate  of 
change  of  "resin"  density,  pf  is  the  time  rate  of  change  of  "fiber"  density, 

Qp  is  the  heat  of  pyrolysis  of  the  resin,  Hy  is  the  heat  of  vaporization  of 
stable  products,  ihg  is  the  mass  flux  of  gaseous  pyrolysis  products,  Cg  is  the 
specific  heat  of  gaseous  pyrolysis  products,  and  I  is  the  laser  irradiance. 

In  applying  this  code  to  tissue,  the  most  straightforward  approach  has  been  to 
treat  the  water  as  a  "fiber"  component  and  any  non-aqueous  components  as 
"resin".  As  mentioned  earlier,  the  term  "fiber"  is  used  here  to  represent  any 
homogeneous  component  of  the  material  which  is  non-decomposing  (but  can 
undergo  a  phase  change).  For  the  present  studies,  however,  tissue  has  been 
treated  as  if  consisting  of  100  percent  "fiber"  having  the  thermophysical 
properties  of  water.  In  future  work,  a  multicomponent  treatment  is  planned  in 
which  the  non-aqueous  constituents  will  also  be  modeled  and  treated  as  a 
"resin"  component  that  can  undergo  pyrolysis  and  charring. 


The  inputs  required  by  this  model  are:  the  target  geometry;  the  target 
boundary  and  initial  conditions;  vaporization  and  pyrolysis  rate  equations  (if 
pyrolysis  is  modeled);  temperature-dependent  thermophysical  properties; 
reflectance;  laser  penetration  depth  in  the  material  as  a  function  of  space 
and  temperature;  and  the  intensity  and  spatial  profile  of  the  laser 
irradiation.  The  effect  of  laser  wavelength  is  treated  by  an  appropriate 
choice  of  tissue  reflectance  and  laser  penetration  depth. 

The  energy  equation  is  solved  with  two  (front  and  back)  boundary  condi¬ 
tions  and  one  initial  condition  using  an  explicit,  forward-marching,  finite- 
difference  technique.  The  (moving)  coordinate  system  is  located  on  the 
receding  front  surface  and  the  numerical  scheme  is  optimized  by  means  of  a 
grid  which  expands  in  time  and  sheds  grid  points  as  the  material  decreases  in 
thickness. 

The  outputs  that  can  be  provided  by  the  model  are:  the  transient  mass 
removal  and  pyrolysis  rates;  the  transient  "resin"  and  char  density;  the 
transient  temperature  distribution;  and  an  energy  balance  of  the  target 
material,  i.e.,  thermal  energy  residing  in  the  target  and  energy  lost  to 
vaporization,  reradiation,  etc. 

2.2.2  Three-Dimensional  Model 


The  three-dimensional  model  accounts  for  transient  heat  conduction  in 
two-  (radial  symmetry)  or  three-dimensional  cartesian  or  cylindrical  polar 
coordinates.  This  model  can  treat  a  variety  of  features,  including  aniso¬ 
tropic  materials,  temperature-dependent  thermophysical  properties,  in-depth 
energy  absorption,  reflection,  and  surface  radiative  and  convective  heat 
transfer  at  the  boundaries  of  the  material.  At  the  present  time,  surface 
ablation  is  simulated  in  the  three-dimensional  code  by  means  of  a  melting 
process  (i.e.,  fixed  vaporization/melt  temperature).  This  is  in  contrast  to 
the  one-dimensional  model  where  the  vaporization  rate,  not  khe  vaporization 
temperature,  is  an  input. 

The  three-dimensional  code  uses  an  implicit  numerical  technique  which  can 
range  from  the  standard  Crank-Nicholson  formulation  to  the  classical  succes- 


2-5 


sive  overrelaxation  implicit  procedure.  The  inputs  required  by  this  model  are 
essentially  the  same  as  those  required  by  the  one-dimensional  model  described 
earlier.  The  outputs  of  this  model  are  ablation  contours  at  various  times  and 
the  transient-temperature  distributions  in  the  target  material. 

2.2.3  Code  Modifications  for  Present  Studies 


In  general,  the  codes  were  used  with  minimal  modifications  with  only  the 
input  changed  to  reflect  the  thermophysical  properties  of  biological  tissue 
(see  discussion  that  follows  on  tissue  material  properties).  The  only  new 
physical  mechanism  added  that  was  not  previously  in  the  codes  is  the  reaction 
that  gives  rise  to  the  thermal  damage  or  necrosis  of  tissue,  i.e.,  denatura- 
tion  of  tissue  protein.  Following  the  work  of  Henriques^  and  others, *^-12  an 
equation  for  damage  was  employed  for  a  cylindrical  coordinate  system  of  the 
form: 


where 


a<r,z)  =  A 


/ 


f 


exp{-E  /RT}  dt 

cl 


(2) 


r,z 

A 

Ea 

R 

T 

Q(r,z) 

*i 

tf 


cylindrical  coordinates 

preexponential  factor  (s-*) 

activation  energy  for  the  reaction  (J/mole) 

universal  gas  constant  =  8.3  J/mole  K 

T(r,z,t)  =  tissue  temperature  (K) 

damage  integral 

onset  time  of  laser  exposure 

time  after  laser  exposure  when  tissue  recovers  to  its  ambient 
temperature 


The  above  damage  equation  is  expressed  in  the  familar  Arrhenius  form  used 
throughout  the  literature  on  thermally-induced  tissue  damage.  The  specific 
values  for  A  and  Ea  used  in  the  calculation  and  the  basis  for  their  selection 
are  presented  in  the  following  subsection. 
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2.2.4  Tissue  Material  Properties 


For  this  initial  modeling  effort,  biological  tissue  has  been  modeled  as  a 
single  component  material  with  constant  thermal  and  optical  properties  similar 
to  those  of  water.  The  philosophy  was  to  begin  our  analysis  with  a  nominal 
baseline  set  of  material  properties  and  then  to  vary  the  laser  irradiation 
conditions  to  explore  the  sensitivity  of  the  results  to  laser  parameters.  In 
general,  the  properties  we  have  used  should  be  representative  of  high  water 
content  tissue,  i.e.,  mucosa, dermal  layers, ^  etc.  The  validity  of  our 
thermal  properties  has  been  tested,  at  least  in  part,  by  comparing  our  model 
predictions  with  a  fairly  extensive  body  of  data  on  laser-induced  thermal 
damage  in  pig  skin.^  The  results  of  those  comparisons  are  presented  in  the 
following  section.  Table  2.1  presents  the  values  of  the  tissue  materials 
properties  used  in  our  calculations. 

TABLE  2.1.  Tissue  Thermophysical  Properties  Used  for  Model  Calculations. 


Thermal/Qptical  Properties 

Density,  p  =  1.05  g/cc 

Specific  Heat,  C  =  4  J/g-K 

Thermal  Conductivity,  k  =  6. 8x10“  8  Vcm-K; 

Thermal  Diffusivity,  K  =  k/pC  =  1.6x10“^  cm^/s 
Heat  of  Vaporization,  Hy  =  2257  J/g 

Absorption  Coefficient  (X  =  10. 6y)  =  769  cm-*  (  a  *  13u) 
Emissivity  (for  thermal  reradiation)  =  1 

Coefficients  in  Denaturation  Equation  (See  Eq.  (2)) 

A  =  3.1xl098  s"1 

Ea  =  150,000  cal/mole  =  628,000  J/mole 

From  Ref.  9 


Future  improvements  in  the  treatment  of  tissue  material  properties  should 
incorporate  material  properties  that  are  a  function  of  the  laser-modified 
water  content. It  is  also  desirable  to  provide  a  realistic  treatment  of 
the  pyrolysis  and  possible  charring  of  the  organic,  non-aqueous  constituents 
of  the  tissue. 

2.2.5  Temperature-Dependent  Vaporization  Rate 

For  the  present  work,  one  improvement  we  have  made  over  previous  analyses 
is  the  incorporation  of  a  tissue  vaporization  rate  that  is  a  function  of 
temperature.  This  has  been  accomplished  by  use  of  a  mathematical  fit  to  the 
equilibrium  vapor  pressure  of  water  for  various  temperatures**  and  coupling 
the  resulting  expression  to  a  detailed  treatment  of  the  Knudsen  layer*  fluid 
dynamics  of  a  surface  vaporizing  into  a  one  atmosphere  background. *^  Thus, 
the  surface  vaporization  temperature  reached  during  a  laser  exposure  is 
actually  a  function  of  the  absorbed  laser  flux,  with  higher  laser  fluxes 
supporting  higher  mass  vaporization  rates,  and  consequently,  higher  surface 
temperatures.  This  improvement  yields  more  realistic  tissue  temperature 
profiles  (surface  temperatures  which  can  exceed  100°C)  and,  as  a  result, 
should  give  more  accurate  predictions  of  the  resulting  subsurface  tissue 
damage . 

2.3  Model  Predictions  and  Comparisons  with  Available  Data 
2.3.1  Low  Power  CV  Exposures 

As  an  initial  evaluation  of  the  basic  thermal  response  model,  we  first 
tested  it  against  available  tissue  response  data  taken  under  low  power  condi¬ 
tions.  Under  these  conditions,  which  are  below  the  threshold  for  significant 


*  The  Knudsen  layer  is  a  narrow  region  just  above  the  vaporizing  surface  (a 
few  mean  free  paths)  in  which  the  gas  particles  undergo  sufficient  collisions 
to  achieve  translational  equilibrium  and  for  the  gas  phase  distribution 
function  to  be  described  by  a  full  Maxwellian. 
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tissue  vaporization,  the  absorbed  laser  flux  goes  into  residual  thermal  energy 
of  the  tissue.  The  resulting  thermal  damage  is  a  function  of  the  tissue 
thermal  properties  only  (density,  specific  heat,  thermal  conductivity),  and 
the  rate  constant  for  the  thermal  denaturation  reaction. 

An  excellent  set  of  data  exists  from  a  series  of  investigations  carried 
out  on  laser-induced  thermal  damage  of  skin.?  In  those  studies,  Takata  and 
covorkers  irradiated  pig  skin  with  CO2  and  ruby  lasers  and  made  histologic 
determinations  of  the  extent  of  thermal  damage.  The  CO2  experiments  were 
performed  with  average  CV  powers  betveen  1.5  and  2V,  Gaussian  beam  profiles, 
and  exposure  durations  from  a  few  tenths  of  a  second  to  several  tens  of 
seconds.  Two  different  beam  spot  sizes  were  employed  with  1/e?  radii  of  0.71 
and  0.383  cm,  respectively,  resulting  in  peak  on  axis  irradiances  of 
1<I0<9  W/cm2. 

In  Figures  2.1  and  2.2,  we  show  typical  comparisons  of  the  predicted 
maximum  depths  of  irreversible  damage  (2  =  0.5+1)12  with  the  corresponding 
observations.7  Our  predictions  are  the  dashed  curves,  the  observations  are 


Figure  2-1.  Model/data  comparisons  of  depths  of  irreversible  tissue 
damage  produced  by  CO2  laser  exposures,  IQ  =  8.55  W/cm2. 
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the  reported  data  points  with  error  bars,  and  the  solid  curves  are  the 
predictions  of  a  model  developed  by  Takata  et  al . ?  As  can  be  seen,  the 
agreement  between  our  model  predictions  and  the  measurements  is  generally 
quite  good,  especially  considering  the  rather  large  error  bars  associated  with 
the  data. 

For  the  calculations  shown,  the  one-dimensional  heat  conduction  code  was 
used  with  the  irradiance  input  being  the  peak  on-axis  value  corresponding  to 
each  experiment  and  the  laser  exposure  duration  also  the  same  as  the  experi¬ 
ment.  Thus,  the  one-dimensional  calculation  was  used  to  predict  the  depth  of 
thermal  damage  occurring  on  axis  (at  the  center  of  the  spot).  Of  course,  the 
measured  maximum  depths  of  damage  also  correspond  to  on-axis  values.  Two- 
dimensional  calculations  run  subsequently  indicated  similar  depths  of  axial 
thermal  damage  clearly  showing  that,  at  least  for  the  range  of  conditions 
modeled,  lateral  heat  conduction  losses  were  not  a  major  factor.  Figure  2.2 
shows  a  comparison  between  the  predictions  of  the  one-  and  two-dimensional 
heat  conduction  codes. 

It  should  be  noted  that  for  most  of  the  laser  exposure  conditions 
modeled,  a  significant  amount  of  thermal  damage  continues  to  occur  even  after 
the  laser  is  shut  off.  This  is  clearly  a  result  of  the  heat  conduction  or 
"thermal  soak"  that  takes  place  subsequent  to  the  laser  exposure.  In  general, 
it  was  found  that  to  assess  the  final  extent  of  thermal  damage  accurately,  it 
was  necessary  to  follow  the  thermal  relaxation  profile  for  at  least  one  full 
exposure  duration  after  termination  of  laser  irradiation. 

The  one-dimensional  calculations  can  predict  fairly  accurately  the  depths 
of  damage  observed  in  the  Takata  experiments  on  axis.  However,  it  is  clear 
that  for  predicting  radial  profiles  of  thermal  damage,  particularly  when 
Gaussian  beams  are  used,  a  two-dimensional  code  is  desirable.  The  results  of 
such  calculations  are  shown  in  Figure  2.3.  As  can  be  seen,  with  a  threshold 
damage  criterion  of  0.1<Qt<l,  the  predicted  and  observed  damage  radii  are  in 
close  agreement. 
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Figure  2-3.  Calculated  radial  profile  of  damage  at  tissue  surface  versus 
experimentally  measured  radius  of  "threshold"  damage. 


2.3.2  Mid  to  High  Power  C V  Exposures 

Calculations  have  been  performed  for  CV  irradiations  ranging  from 
250  V/cm2  to  5  KV/cm2  and  with  exposure  durations  as  short  as  10-2s  and  as 
long  as  10s.  Outputs  of  the  code  included  mass  loss,  surface  recession  (i.e., 
crater  depth),  residual  thermal  energy,  spatial/temporal  temperature  profiles, 
and  profiles  of  the  damage  integral  (Q).  At  least  for  the  higher  irradiances 
considered,  ve  have  found  that  steady  state  conditions  are  predicted  to  occur 
early  in  the  laser  exposure  (within  approximately  3  ms),  at  which  time  the 
surface  temperature,  recession  rate,  and  residual  thermal  energy  all  assume 
constant  values  until  the  laser  shuts  off. 

The  computed  in-depth  temperature  profiles  and  damage  integral  (Q)  con¬ 
tours  are  plotted  in  Figures  2.4  and  2.5  for  an  incident  irradiance  of 
5  ktf/cm2  before  the  laser  shuts  off  and  0.05s  later.  (Note  that  the  depths 

are  given  in  reference  to  the  original  tissue  surface  and  that  there  has  been 

approximately  0.9  mm  of  tissue  recession.)  An  interesting  feature  of  Figure 
2.5  is  that,  while  the  laser  is  on,  the  peak  temperature  is  predicted  to  occur 
beneath  the  surface.  Other  researchers  modeling  laser-induced  vaporization  of 
solid  materials  with  in-depth  absorption  have  predicted  similar  temperature 
profiles. For  example,  in  Ref.  16  it  was  found  that  "for  certain  laser  and 
material  parameters,  subsurface  temperatures  will  significantly  exceed  the 
surface  temperature,  and,  under  such  conditions,  explosive  removal  of  material 
could  occur,  resulting  in  very  rapid  and  efficient  material  removal." 

It  is  interesting  to  note  (from  Figure  2.4)  that  while  there  is 

<1  percent  additional  contribution  to  tissue  recession  after  the  laser  shuts 
off,  the  zone  of  in-depth  thermal  damage  grows  by  about  60  percent  (from 
-50  ym  to  -80  ym).  This  demonstrates  that  modeling  the  in-depth  conduction  of 
residual  heat  after  the  laser  shuts  off  is  important  for  accurately  predicting 
the  extent  of  thermal  damage.  This  will  become  even  more  critical  when 
dealing  with  short  pulse  interactions. 
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Figure  2-4.  Computed  in-depth  damage  profiles  after  CW  irradiation 
of  5  kV/cm^  (laser  exposure  time  is  0.05s). 
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Figure  2-5.  Computed  in-depth  temperature  profiles  in  tissue 
(I0  *  5  kV/cm^,  exposure  duration  *  0.05s). 


In  Table  2.2  we  summarize  the  results  of  four  CV  calculations  run  for 
irradiances  of  250,  1300,  2000,  and  5000  V/cm^  vith  exposure  times  of  1,  0.19, 
0.125,  and  0.05s,  respectively.  The  objective  vas  to  generate  the  same 
exposure  fluence  (250  J/cm^)  at  several  different  irradiance  levels  and  to 
compare  their  relative  effect  on  recession  and  the  extent  of  residual  thermal 
damage.  One  of  the  more  interesting  quantities  in  Table  2.2  is  Zr/ AZj  or  the 
ratio  of  the  recession  depth  to  the  residual  depth  of  irreversible  tissue 
damage.  Clearly,  if  the  goal  is  to  optimize  tissue  cutting  while  minimizing 
the  extent  of  residual  tissue  damage,  then  larger  values  of  Zr/AZ<j  are  desir¬ 
able.  It  is  clear  from  Table  2.2  that  the  smallest  relative  amount  of 
residual  damage  is  achieved  when  operating  at  the  highest  irradiance  level  and 
shortest  exposure  time.  These  results  provide  quantitative  support  for  a 
prevailing  opinion  that  has  now  become  increasingly  evident  in  the  laser 
treatment  literature. 17-19 
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TABLE  2.2.  Comparison  of  CV  Laser  Exposures  at  Different  Irradiance 
Levels  But  Fixed  Exposure  Fluence  (250  J/cm2). 


Exposure 

mm 

■■■ 

Irradiance 

Duration 

AZd+ 

(V/cm2) 

(s) 

(cm) 

Zr/AZd 

250 

1.00 

4.83 

0.1048 

0.0244 

4.3 

1300 

0.192 

2.77 

0.1144 

0.0144 

7.9 

2000 

0.125 

2.44 

0.1133 

0.0125 

9.1 

5000 

0.050 

1.91 

0.0928 

0.009 

10.3 

♦Qr  *  final  tissue  residual  energy 

+Zr  a  final  tissue  recession  depth 

♦AZd  -  final  residual  depth  of  irreversible  tissue 

damage 

It  is  worth  mentioning,  however,  that  there  are  generally  other  factors 
which  set  upper  limits  to  the  irradiance  level  (and,  correspondingly,  lower 
limits  to  the  pulse  duration)  which  is  most  useful  for  clean,  efficient, 
well-controlled  ablation.  One  of  these  factors  is  the  laser  penetration  1 

depth.  For  example,  it  is  expected  that  little  value  will  be  gained  by 
decreasing  the  exposure  duration  to  much  less  than  $2/K  where  &  is  .he  optical  | 

penetration  depth  and  K  is  the  thermal  diffusivity  of  the  material.  Beyond 
this  point,  the  deposition  of  heat  during  the  laser  exposure  time  will  be  | 

controlled  more  by  the  scale  length  of  the  optical  deposition  than  by  thermal 
conduction.  Other  factors  which  may  come  into  play  at  high  irradiance  levels  \ 

are  explosive  material  blow-off  phenomena  and  the  effects  of  plasma  formation. 
Depending  on  the  situation,  such  processes  could  be  undesirable.  . 
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It  should  also  be  mentioned  that  the  assertion  that  operating  at  higher 
irradiances  and  shorter  exposure  times  produces  less  residual  damage  is 
necessarily  true  only  if  the  spot  size  is  kept  constant.  If,  on  the  other 
hand,  the  higher  irradiance  is  achieved  by  focussing  to  a  smaller  diameter 
spot,  radial  heat  conduction  losses  and  the  relative  extent  of  lateral  thermal 
damage  may  increase. 

Figure  2.6  demonstrates  the  ability  of  the  two-dimensional  code  to 
predict  the  growth  of  residual  thermal  damage  around  the  perimeter  of  a  small 
spot  crater.  A  two-dimensional  calculation  was  run  for  a  Gaussian  beam  with  an 
intensity  profile  of  I(r)  =  1990  exp[-(r/0.031)2]  V/cm^  (1/e  radius  =  310  urn) 
and  exposure  duration  of  0.5s.  Figure  2.6  presents  the  computed  results  for 
the  crater  profile,  the  zone  of  residual  thermal  damage,  and  temperatures  at 
selected  points  in-depth  in  the  tissue. 


It  would  be  extremely  valuable  to  test  model  predictions  of  the  type 
shown  in  Figure  2.6  against  comparable  experimental  data.  Unfortunately,  we 
have  not  yet  found  an  appropriate  body  of  data  taken  at  these  higher  fluxes 


Figure  2-6.  Calculated  two-dimensional  crater  profile  and  temperature 
map  at  0.5s  (P0  =  6W,  I  =  1990  exp[ l-(r/0. 031)2 J ) . 
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which  is  of  acceptable  quality.  In  particular,  most  of  the  tissue  cutting 
data  available  in  the  literature  lacks  an  adequate  specification  of  the  laser 
beam  spatial  intensity  profile  that  was  used.  A  major  goal  of  future  efforts 
should  be  to  acquire  well-characterized  and  well-diagnosed  laser/tissue 
ablation  data  to  serve  as  a  rigorous  test  of  the  models  that  have  been 
developed  and  to  determine  where  improvements  are  needed. 

2.3.3  Modeling  of  Short,  Single-Pulse  Exposures 

It  is  possible  with  present  surgical  CO2  lasers  to  deliver  laser  energy 
to  a  treatment  site  in  a  variety  of  enhanced  pulsing  modes.  In  these  modes, 
the  laser  is  electronically  pulsed  to  generate  higher  peak  powers  (up  to  a 
factor  of  5  to  7  above  the  average  CV  power)  in  short  repetitive  bursts 
(presently  10“^  to  10-^s  in  duration).  A  considerable  body  of  data  in  the 
laser  surgical  literature  indicates  that  this  enhanced  or  "super-pulsed"  mode 
of  operation  can  yield  cleaner  cuts  with  less  residual  damage  to  adjacent 
tissue.  It  should  be  mentioned,  however,  that  these  assertions  are  not 
without  controversy.  In  addition,  there  is  now  a  developing  data  base  for 
which  pulsed  operation  is  being  extended  down  to  the  microsecond  (10~6s) 
regime. 8»19  Based  on  preliminary  analyses  of  the  latter  experiments,  even 
more  interest  is  being  generated  in  the  potential  of  using  short  pulses  for 
achieving  "clean"  tissue  cutting. 

Code  Predictions 


To  address  some  of  these  issues,  the  one-dimensional  ablation  code  was 
used  to  run  a  series  of  single-pulse  calculations  for  a  range  of  pulse 
fluences  varying  from  1  to  100  J/cm^  per  pulse  and  for  pulse  durations  from 
10-5  to  10“^s.  Typical  results  are  presented  in  Figures  2.7  and  2.8. 

Figure  2.7  is  a  plot  of  the  calculated  residual  thermal  energy  in  the 
tissue  just  after  the  ablation  stops.  The  ratio  of  this  quantity  to  the 
incident  pulse  fluence  is  the  thermal  coupling  coefficient.  Note  that  at  the 
higher  pulse  fluences,  i.e.,  =  100  J/cm^,  the  thermal  coupling  coefficient  is 
predicted  to  be  small,  -1  to  2  percent.  The  reason  is  that  the  energy 
deposition  rate  is  so  high  that  steady  state  vaporization  of  the  tissue  is 


Figure  2-7.  Calculated  residual  thermal  energy  in  tissue  as  a  function 
of  pulse  fluence  for  several  different  pulse  durations. 


achieved  early  in  the  pulse,  so  that  the  remainder  of  the  pulse  energy  is 
simply  carried  off  by  vaporization  products. 


In  Figure  2.8  we  present  the  corresponding  calculations  of  the  recession 
depth  and  depth  of  thermal  damage  (relative  to  the  original  surface  location). 
The  depth  of  residual  tissue  damage  is  simply  the  difference  between  the 
latter  and  the  former.  It  is  clear  from  this  plot  that  the  ratio  of  the  depth 
of  residual  tissue  damage  to  the  recession  depth  is  decreased  by  operating  at 
the  higher  fluences. 


The  results  of  these  preliminary  single-pulse  calculations  have  led  us  to 
the  following  conclusions: 

The  most  effective  tissue  cutting  with  minimal  residual  thermal  energy 
and,  hence,  minimal  residual  thermal  damage,  appears  to  be  best  achieved 
by  operating  with  pulses  i  10'^s  in  duration  and  at  pulse  fluences 
>  30  J/cm^,  or,  with  peak  pulse  ir radiances  >  300  kW/cm^.  Under  these 
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DEPTH  MEASURED  FROM  ORIGINAL  SURFACE  (cm) 


Figure  2-8.  Calculated  tissue  recession  depth  and  thermal  damage  depth 
(relative  to  original  surface)  as  a  function  of  pulse 
fluence  for  several  different  pulse  durations. 


conditions,  tissue  recession  velocities  are  predicted  to  be  so  rapid  that  the 
recession  front  literally  keeps  up  with  the  thermal  wave  front  and,  as  a 
result,  inhibits  the  deposition  of  heat  beneath  the  receding  surface. 

Analytical  Blow-Off  Model 

One  simple  thermal  model  of  pulsed  laser  ablation  which  is  occasionally 
employed  is  a  critical  energy  "blow-off"  model  based  on  Beer's  law.  This 
model  assumes  the  laser  energy  fluence  is  deposited  in  depth  according  to 
Beer's  law,  i.e.,  F(Z)  -  F0e-2/$  (where  S  is  the  (1/e)  absorption  depth  at  the 
laser  wavelength),  and  that  the  target  material  is  ablated  to  a  depth  at  which 
the  deposited  energy  just  falls  below  a  critical  energy  density  required  for 
ablation.  Because  it  neglects  thermal  conduction,  this  model  can  only  be 
applied  to  short  pulse  exposures  for  which  tp  <  S^/K,  where  k  is  the  thermal 
diffusivity  of  the  target  material.  For  example,  for  tissue  interactions  at 
CO2  wavelengths  (A  =  10.6  microns,  S  =  1.3  x  10“^  cm),  this  would  correspond 
to  pulse  durations  <  10“  ^  seconds.  The  above  assumptions  lead  to  the  follow¬ 
ing  relationships: 


9F 

9Z 

Z=Z  =  ec  -  (1-R)  V5  e  r 

(3) 

r 

fT  = 

Sec/1-R 

(M 

IS 

II 

6  n(F0/FT) 

(5) 

and 


where 


^residual  -  FT 


ec  =  critical  energy  density  for  ablation  (J/cc) 
=  P  habl  =  P(hv  +  C(TV-T0)) 

Zr  =  material  recession  or  "blow-off"  depth 
R  =  target  reflectance 

Fj  =  threshold  fluence  for  ablation  (J/cm^) 
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^residual  =  thermal  energy  remaining  in  target  after  ablation  ends 
(J/cm2). 

For  tissue  (water),  hv  is  taken  to  be  2257  J/g,  C  =  4  J/g°C,  Tv  =  100°C,  and 
T0  =  37°C.  This  gives  ha5j  *  2.6  kJ/g. 

Belov,  we  apply  this  model  to  the  pulsed  tissue  ablation  experiments  of 
Valsh  et  al.  and  compare  the  results  to  corresponding  predictions  of  the 
numerical  heat  conduction  model. 

Comparisons  with  Data 

The  only  tissue  ablation  data  we  found  that  could  reasonably  be  used  to 
assess  single-pulse  tissue  response  is  the  data  of  Walsh  et  al.,®  recently 
obtained  at  the  Wellman  Laboratory,  Massachusetts  General  Hospital.  In  those 
experiments,  the  ablation  and  residual  damage  of  pig  skin  was  studied  using 
2-us  long  pulses  from  a  CO2  TEA  laser.  A  1.6  x  2.0  mm  spot  size  was  used  in 
all  experiments.  The  residual  tissue  damage  was  assessed  from  laser  exposures 
performed  in  vivo,  after  which  the  irradiated  site  was  biopsied  and  the  tissue 
processed  for  histological  examination.  Mass  removal  rates,  on  the  other 
hand,  were  measured  in  vitro  by  mounting  a  sample  of  pig  skin  on  an  electronic 
balance  and  recording  the  weight  loss  resulting  from  laser  exposure  using  a 
computer-controlled  data  acquisition  system. 

Walsh's  observations  of  mass  loss  are  presented  in  Figure  2.9.  The 
single-pulse  mass  loss  is  inferred  in  these  experiments  by  accumulating  a 
train  of  pulses  on  one  tissue  sample  at  a  low  repetition  rate  (approximately 
2  Hz),  measuring  the  accumulated  mass  loss,  and  effectively  dividing  by  the 
number  of  pulses.  The  number  of  pulses  delivered  varied  from  40  at  high 
fluences  to  130  at  the  lowest  fluence.  Plotted  in  Figure  2.9  is  the  average 
mass  ablated  per  pulse  per  unit  area  versus  the  average  incident  fluence  per 
pulse  (J/cm2/pulse). 

Also  plotted  with  the  data  are  corresponding  single-pulse  predictions  of 
the  one-dimensional  heat  conduction  code  (solid  curves)  and  the  predictions  of 
the  Beer's  law  blow-off  model  (dashed  curves)  for  two  different  assumed  values 
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Figure  2-9.  Model/data  comparisons  of  tissue  mass  loss  per  pulse 
versus  fluence  per  pulse  for  irradiation  by  micro¬ 
second  duration  CO2  TEA  laser  pulses. 

of  the  absorption  depth.  The  (1/e)  absorption  depth  previously  used  in  the 
thermal  conduction  calculations  was  13  microns.  The  calculated  areal  mass 
loss  is  related  to  recession  depth  simply  as  dim  =  pZr  where  p  is  the  tissue 
density  (1.05  g/cc  used). 

There  are  several  significant  features  of  the  results  shown  in  Figure  2.9 
that  are  worth  noting.  Perhaps  the  most  important  feature  is  that  the 
observed  dependence  of  the  mass  loss  on  pulse  fluence  appears  to  be  more 
linear  than  logarithmic  and,  therefore,  fails  to  follow  the  predictions  of  a 
blow-off  model.  Even  if  a  shorter  absorption  depth  or  lower  critical  energy 
for  vaporization  were  arbitrarily  assumed  to  better  fit  the  data  at  the  lower 
fluences,  there  would  then  be  an  even  larger  discrepancy  at  the  higher 
fluences  (the  blow-off  model  would  significantly  underpredict  the  mass  loss). 

The  dependence  of  the  measured  mass  loss  on  pulse  fluence  is  clearly  in 
better  agreement  with  that  predicted  by  the  thermal  conduction  model.  In 
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addition,  the  thermal  conduction  model  much  more  closely  predicts  the  fluence 
"threshold"  at  which  measurable  ablation  is  observed  (-1  J/cm2).  The  only 
significant  discrepancy  between  the  predictions  of  the  thermal  conduction 
calculations  and  the  data  is  that  the  measured  mass  loss  is  generally  about  20 
to  25  percent  lower  than  the  predictions.  We  offer  a  possible  explanation  for 
this  discrepancy  at  the  conclusion  of  this  discussion. 

This  apparent  linear  dependence  of  mass  loss  on  pulse  fluence  suggests 
that,  at  least  for  pulse  fluences  sufficiently  above  "threshold,"  a  model  of 
"steady  state”  vaporization  is  much  more  appropriate  than  a  critical  energy 
blow-off  model.  It  is  also  significant  that  measurable  mass  loss  is  observed 
at  a  lower  fluence  than  the  critical  threshold  fluence  defined  in  the  blow-off 
model.  This  is  consistent  with  a  thermal  conduction  and  surface  vaporization 
model  which  allows  vaporization  to  occur  as  soon  as  the  front  surface 
temperature  reaches  the  boiling  point,  i.e.,  Fabs  =  pC(Tb~T0).&  =  0.34  J/cm2 
(for  5  =  13  urn).  Subsequent  to  that  point,  the  mass  vaporization  rate  (and 
corresponding  surface  temperature)  during  the  pulse  are  determined  by  the 
absorbed  flux  rather  than  explicitly  by  fluence. 

In  a  similar  manner,  the  assumptions  of  the  blow-off  model  lead  to  an 
overestimate  of  the  residual  thermal  energy  or  retained  heat  in  the  tissue 
after  ablation  terminates.  For  example,  we  see  from  Figure  2.7  that  for  an 
incident  fluence  of  10  J/cm2,  the  thermal  conduction  model  predicts  a  final 
residual  thermal  energy  of  -  1.2  J/cm2  while,  for  the  same  absorption  depth  of 
13  pm,  the  blow-off  model  predicts  -  3.3  J/cm2,  or  about  a  factor  of  3  higher. 
Part  of  the  reason  for  this  difference  is  that  in  the  conduction  model  some 
evaporation  can  continue  for  a  significant  time  after  the  laser  pulse 
terminates,  allowing  the  tissue  to  cool  evaporatively  and  leaving  less  heat 
behind,  while  in  the  blow-off  model  no  evaporation  is  permitted  after  the 
pulse  terminates.  In  the  thermal  conduction/surface  vaporization  model, 
ablation  will  continue  until  the  front  surface  temperature  drops  to  ~  100°C 
(or  local  incremental  energy  density  =  260  J/g).  This  is  nearly  one-tenth  of 
the  critical  energy  density  at  which  vaporization  is  forced  to  terminate  in 
the  blow-off  model. 
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These  differences  in  the  residual  thermal  energies  predicted  by  the  tvo 
models  also  lead  to  differences  in  the  extent  of  residual  thermal  tissue 
damage  that  is  predicted  to  occur.  To  first  order,  in  fact,  one  vould  expect 
the  zone  of  residual  damage  to  scale  proportionally  vith  the  residual  thermal 
energy.  For  example,  at  an  incident  fluence  of  10  J/cm^  and  S  -  13  pm,  we 
calculate  from  the  code  Qr  =  1.2  J/cm^  and  a  final  depth  of  residual  damage  of 
51  microns.  In  contrast,  if  ve  allow  the  corresonding  higher  residual  energy 
of  the  blow-off  model  (Qr  =3.3  J/cm^)  to  conduct  into  the  tissue  and  calcu¬ 
late  the  denaturation  which  results,  we  find  final  residual  damage  depths  of 
between  130  to  160  microns. 

The  observed  depths  of  residual  damage  in  the  Valsh  experiments  were 
quoted  to  be  between  50  and  100  microns  over  the  pulse  fluence  range  of  2  to 
18  J/cm^.  For  the  code  calculations  carried  out  over  the  same  fluence  range 
(and  for  an  assumed  absorption  depth  of  between  13  and  20  microns),  damage 
depths  are  found  to  be  between  40  and  80  microns.  On  the  other  hand,  using 
the  corresponding  residual  energies  given  by  the  blow-off  model,  damage  depths 
are  predicted  to  be  from  a  minimum  of  130  microns  to  a  maximum  of  240  microns. 
Thus,  once  again,  we  find  the  code  predictions  to  be  in  closer  agreement  with 
the  observations. 

Finally,  we  have  considered  why  the  measured  "single-pulse"  mass  loss  in 
Valsh's  experiments  is  generally  lower  than  the  single-pulse  code  predictions. 
At  this  point  we  believe  the  most  likely  explanation  comes  from  the  fact  that 
the  experimental  "single-pulse"  mass  loss  was  obtained  by  averaging  over  the 
mass  loss  from  an  accumulation  of  multiple  single-pulse  exposures  to  the  same 
area  of  tissue.  It  is  therefore  quite  likely  that,  in  the  process,  the 
surface  of  the  tissue  became  dried  out  or  dessicated,  thus  leaving  less  tissue 
water  available  for  vaporization  on  subsequent  pulses.  Ve  suggest  that  in  the 
future  an  attempt  be  made  to  obtain  true  "single-pulse”  data.  In  addition, 
modeling  efforts  should  be  directed  toward  a  multicomponent  treatment  of 
tissue,  i.e.,  aqueous  and  non-aqueous  components.  This  would  allow  for  an 
accounting  of  the  depletion  of  available  tissue  water.  A  proper  treatment, 
however,  will  require  a  model  for  the  pyrolysis  and  possible  charring  of  the 
non-aqueous  residue. 
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2.4  Summary  and  Conclusions 


Detailed  model  calculations  have  been  carried  out  to  assess  the  thermal 
response  and  ablation  of  biological  tissue  exposed  to  continuous  wave  and 
short  pulse  CO2  laser  radiation.  The  models  developed,  however,  are 
applicable  to  any  laser  wavelength  for  which  absorption  in  tissue  dominates 
over  scattering. 

The  model  predictions  are  shown  to  be  in  close  agreement  with  available 
data  on  tissue  damage  and  ablation.  The  results  confirm  that  a  transient 
thermal  conduction  model  that  includes  an  equation  for  the  thermal  denatur- 
ation  of  tissue  protein  and  the  surface  vaporization  of  tissue  water  can  in 
large  measure  explain  the  experimental  observations. 

Future  modeling  improvements  that  are  planned  are  a  multicomponent 
treatment  of  tissue  that  will  model  the  non-aqueous  as  well  as  the  aqueous 
components,  and  the  inclusion  of  a  laser  source  term  that  accounts  for 
scattering. 
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3.  PLASMA  MODELING 


Experiments  on  laser-induced  fragmentation  of  biliary  and  urinary  calculi 
have  shown  that  fragmentation  is  alvays  accompanied  by  plasma  formation^ * 2 . 

The  plasma  is  visible  as  a  bright  flash  of  duration  of  the  order  of  the  pulse 
length.  Spectroscopic  investigation  of  the  radiation  emitted  by  biliary 
calculi  during  the  flash  has  shovn  that  the  emission  consists  of  a  continuum 
spanning  the  visible  spectrum  with  strong  lines  of  ionized  Calcium  appearing 
in  absorption2.  The  continuum  disappears  at  a  later  time  and  the  lines  then 
appear  as  emission  features. 

Ve  have  investigated  mechanisms  consistent  vith  the  experimental  data 
through  which  a  plasma  could  be  ignited  in  the  calculi.  The  data  indicate 
that  the  threshold  is  determined  primarily  by  the  laser  fluence  rather  than  by 
the  laser  irradiance.  Graham  Vatson1  performed  experiments  using  various 
lasers  with  pulse  lengths  extending  from  20  ns  (Q  switched  ruby  and  Nd  Yag 
lasers)  to  300  ys  (flash  lamp  pumped  dye  laser).  He  quotes  at  the  shortest 
pulse  length  a  threshhold  of  40  mJ  at  the  doubled  Yag  frequency  (X  *  0.53  ym) 
with  a  beam  diameter  at  focus  of  0.5  mm.  This  corresponds  to  a  fluence  4  » 

40  x  10-3  /  (n  x  (0.05)2  /  4)  «  20  J/cm2.  This  number  is  to  be  compared  with 
a  fluence  of  34  J/cm2  at  the  longest  pulse  length  studied,  t  »  300  ys, 
corresponding  to  the  conditions  X  ■  0.577  ym,  d  *  1  mm.  The  threshold 
irradiances,  however,  varied  over  several  orders  of  magnitude  (I  -  10^  V/cm2 
for  t  «  20  ns  and  I  »  10^  V/cm2  for  x  -  300  ys).  The  fluence  thresholds, 
which  are  found  to  lie  in  the  range  4-40  J/cm2,  appear  also  to  be  relatively 
insensitive  to  laser  wavelength. 

Figure  3.1  shows  plasma  flash  data  taken  by  Teng  et  al2  who  used  a  dye 
laser  (X  »  690  nm,  x  -  0.8  ys)  on  a  pigmented  biliary  calculus.  This 
figure  shows  that  breakdown  has  occurred  within  400  ns  after  the  beginning  of 
the  pulse  and  that,  at  the  onset  of  breakdown,  5  J/cm2  have  been  deposited 
in  the  stone.  The  beam  irradiance  is  estimated  to  be  I  *  8  J/cm2  /  600  ns  * 
107  W/cra2. 

Laser-induced  breakdown  in  semi-transparent  optical  materials  has  been 
studied  in  great  detail^.  The  threshold  irradiances  in  pure  crystalline 
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Figure  3.1.  Time  dependence  of  laser  and  flash  intensities. 

Arrovs  mark  the  beginning  of  the  laser  pulse 
and  of  plasma  emission  (from  Teng  et  al*). 

materials  devoid  of  impurities  and  defects  are  expected  to  be  quite  high,  in 
excess  of  101®  W/cm^.  These  thresholds  are  much  higher  than  those  that  have 
been  observed  for  the  breakup  of  stones.  Much  lover  breakdovn  thresholds  have 
been  observed  in  optical  materials  due  to  surface  or  bulk  defects.  The 
surface  defects  can  be  cracks,  grooves,  or  spherical  pores4.  There  may  also 
be  dust  particles  or  surface  inclusions.  Bulk  impurities  could  be  platinum 
particles  from  the  crucibles  used  in  the  fabrication  process.  In  a  similar 
fashion,  absorbing  inclusions  in  biliary  calculi  would  come  from  the  presence 
of  pigment.  Absorption  sites  will  result  in  localized  heating  that  under 
certain  conditions  can  lead  to  thermal  runaway. 

It  seems  apparent,  from  the  low  irradiance  thresholds  observed  for  break¬ 
dovn  on  gall  stones  and  urinary  tract  calculi,  that  plasma  ignition  is  associ¬ 
ated  vith  absorption  and  ignition  on  localized  sites.  The  stones  have  long 
absorption  lengths  in  the  visible  which  means  that  the  radiation  can  penetrate 
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deeply.  At  the  breakdown  fluences  of  tens  of  joules  per  square  centimeter, 
one  does  not  expect  sufficient  bulk  heating  to  account  for  the  phenomena 
observed.  Vhat  ve  believe  is  happening  is  that  the  temperature  around  the 
absorption  site  is  high  enough  to  trigger  non-linear  absorption  in  the 
surrounding  medium.  If  this  absorption  is  high  enough,  the  added  heating  will 
lead  to  thermal  runaway.  There  will  be  a  threshold  irradiance  associated  with 
the  minimum  temperature  to  initiate  the  runaway  and  a  threshold  fluence  asso¬ 
ciated  with  the  development  of  individual  plasmas  at  each  site  and  coalescence 
of  the  plasmas.  Plasmas  will  radiate  in  the  UV  and  deep  UV.  Absorption  of 
these  high-energy  photons  by  the  calculus  will  cause  ionization  of  the 
surrounding  medium.  Inverse  Bremsstrahlung  absorption  by  the  electrons  thus 
generated  will  result  in  enhanced  absorption  and  provides  a  mechanism  for 
plasma  growth. 

Absorption  by  inclusions  of  various  sizes  has  been  analyzed  by  Sparks  and 
Outhler^.  These  authors  considered  temperature  rises  of  only  1000  K,  since 
they  were  solely  concerned  with  material  failure  due  to  stress  buildup. 
Anisimov  and  Makshantsev  also  considered  absorption  by  inclusions  and  postu¬ 
lated  a  thermally-dependent  absorption  coefficient  of  the  medium  of  the  form6: 

k  =  b  exp  -  (E/T)  (1) 

where  E  is  an  energy  of  the  order  of  the  band  gap  of  the  material.  They  were 
able  to  calculate  the  threshold  irradiance  as  a  function  of  the  particle  size. 
We  review  their  theory  below. 

Theory  of  Anisimov  and  Makshantsev. 

Let  there  be  an  absorbing  inclusion  of  radius  R,  embedded  in  a  non¬ 
absorbing  medium.  After  a  transient  period,  a  steady  state  temperature 
profile  will  be  established  in  which  the  flow  of  heat  out  of  the  inclusion 
balances  the  power  absorbed  from  the  laser  field.  If  we  assume  a  spherically 
symmetric  temperature  profile,  then  the  heat  conduction  equation  in  the 
surrounding  medium  can  be  written  as 
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subject  to  the  boundary  condition  T  *  T®  when  r  -»  ®  and,  at  r  =  R, 
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where  I  is  the  laser  irradiance,  K  the  thermal  conductivity,  and  Qabs  the 
absorption  efficiency  of  the  inclusion.  The  solution  to  Eqs.  (2)  and  (3), 
assuming  that  K  is  independent  of  T,  is: 


T  -  T 


Qabs  I  R2 


Anisimov  and  Makshantsev  assume  that  the  medium  surrounding  the  inclusion 
has  a  temperature-dependent  absorption  coefficient  of  the  form  given  by 
Eq.  (1).  Such  a  form  is  reasonable  if  one  assumes  that  a)  free  or  bound  elec¬ 
tronic  states  can  be  thermally  populated  into  allowed  states  having  energy  E 
above  the  Fermi  level  and  b)  the  electrons  in  these  states  can  absorb  the 
laser  photons.  Because  of  this  extra  absorption,  the  temperature  will  rise 
above  the  solution  given  by  Eq.  (4).  In  order  to  calculate  this  temperature 
rise,  one  must  solve  the  time-dependent  heat  equation: 


3T  1  3  r  2  3T^ 

^ St - 2  lr  K2Fj 


where  p  and  C  are  the  density  and  specific  heat  of  the  material.  A  solution 
to  (5)  can  be  obtained  by  expanding  around  the  time-independent  solution  T0(r) 
which  is  approximated  by  Eq.  (4)  (since  k  is  assumed  small) 

T(t)  -  T0(r)  +  0(t,r) 


and  writing  that 


3-4 


k(T) 

We  expand  6  as  follows 
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Inserting  the  expansion  (6)  into  Eq.  (5),  one  finds  that  each  Yy  should  satisfy 

(7) 


4-  -t?  <r2v  +  VY  +  V(r))*  0 

r 


where 


and 


pCy 

Y  "  TT 


V(r)  « 


Anisimov  and  Hakshantsev  seek  a  solution  to  Eq.  (7)  with  the  boundary 
conditions 

jy 

Y  -»  0  when  r  -»  ®  and  -g— 

the  latter  condition  meaning  that  the  additional  heat  absorption  by  the  medium 
surrounding  the  inclusion  does  not  affect  the  heat  flow  out  of  the  inclusion 
(which  is  just  the  power  absorbed  by  the  inclusion  from  the  laser  beam).  If  a 
solution  to  Eq.  (7)  can  be  found,  that  admits  a  negative  eigenvalue  y  (or  y), 
then,  following  Eq.  (6),  one  sees  that  there  will  be  a  thermal  runaway.  Using 
the  fact  that  V(r)  is  appreciably  different  from  zero  only  in  the  immediate 
vicinity  of  R,  they  find,  by  integrating  Eq.  (7)  over  r  once  and  using  the 
boundary  condition  on  Y,  that  a  solution  with  negative  y  can  exist  only  if  the 
following  condition  on  k  is  satisfied. 


J  k'  dr  > 


K 

TR 


(8) 


Now, 


a 

i 


k'dr 


T(R) 


rT»  dk 

(Jr  ) 

J  dT 

(  dT  J 

dT  « 


k[To(R)] 


(dT/dr) 


-1 

r=R 
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If  we  now  use  Eq.  (3)  to  calculate  K(dT/dr)R,  we  can  recast  Inequality  (8)  as 
follows: 


k(T0)  > 


^abs 


or,  using  Eq.  (4)  for  (T0  -  Tod), 


(9) 


Equation  (9)  is  the  fundamental  result  of  Anisimov  and  Makshantsev.  It 
differs  slightly  from  their  Eq.  (8)  in  which  they  dropped  the  T®  term  and 
where,  due  to  a  misprint,  the  logarithm  term  appeared  in  the  numerator  instead 
of  the  denominator.  Also,  rather  than  using  Qabs,  they  defined  an  absorption 
factor  a  (*  Qabs/4)* 


Discussion 


Inequality  (9)  shows  the  breakdown  threshold  irradiance  scales  roughly  as 
R-l.  There  will  be  a  whole  range  of  sizes  for  the  defects  in  the  focal  volume 
and  the  larger  ones  will  be  the  most  easy  to  break  down.  In  a  given  experi¬ 
ment,  the  irradiance  I  is  specified;  therefore,  one  can  expect  that  all 
inclusions  whose  radius  exceeds  a  threshold  value  corresponding  to  Eq.  (9), 
being  satisfied,  will  be  a  plasma  ignition  site.  In  order  to  test  whether  the 
theory  of  Anisimov  and  Makshantsev  is  consistent  with  the  experimental  data, 
we  must  a)  show  that  Inequality  (9)  allows  low  threshold  irradiance  and 
b)  explain  why  there  is  also  an  important  effect  of  fluence.  As  a  typical 
example,  let  us  set  E  =  60,000  K  corresponding  to  a  band  gap  of  5.5  eV,  Qabs  = 
0.5,  K  ■  10~2  V/cmK  and  a  limiting  high  temperature  absorption  coefficient  k  = 
105  cm-*  (i.e.,  b  *  10^).  We  find  that  for  R  =  1  Mm,  I  =  2  x  107  U/cm^, 

while  for  R  =  10  ym,  I  *  7  x  10^  W/cm^.  These  thresholds  are  not  unreasonable. 

The  fact  that  there  is  also  a  threshold  fluence  that  is  weakly  dependent 
on  pulse  length  can  arise  from  two  phenomena.  First,  one  must  deposit  enough 
energy  in  the  inclusion  for  the  steady  state  profile  to  be  established. 

Second,  the  plasmas,  once  ignited,  require  a  certain  time  to  grow  and  coalesce. 
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Estimate  of  Laser  Ignition  Time  and  Threshold  Fluence 


Ve  can  obtain  a  rough  estimate  of  the  heating  time  by  looking  at  the 
energy  balance  equation.  The  heating  curve  is  shown  in  Figure  3.2.  There  is 
an  initial  transient  that  lasts  a  time  Tg  during  which  the  temperature  of  the 
inclusion  rises  linearly  with  time.  After  this  time,  the  rate  of  temperature 
rise  decreases,  the  temperature  asymptoting  to  a  steady  state  value  Ts.  The 
departure  from  the  linear  behavior  occurs  when  the  thermal  diffusion  length 
into  the  surrounding  medium,  Iq,  becomes  comparable  to  the  radius  of  the 
inclusion 


/D  =  /xtH  =  0.2  to  0.5  x  R 


where  X  is  the  thermal  diffusivity  and  tfl  the  thermal  diffusion  time.  The 
thermal  runaway  that  was  discussed  in  the  previous  paragraphs  can  occur  under 
two  different  regimes,  depending  on  the  irradiance  level.  If  ytn  »  1,  then 
thermal  runaway  occurs  during  the  rising  portion  of  the  temperature  versus 


Figure  3.2.  Temperature  versus  time  at  the  surface  of 
an  absorbing  inclusion. 


time  curve,  i.e.,  breakdovn  is  initiated  before  the  particle  temperature  can 
attain  a  steady  state.  In  the  opposite  limit,  rtf]  <  1,  thermal  runaway  will 
occur  long  after  steady  state  heating  has  been  achieved.  Due  to  the  strong 
temperature  dependence  of  k  on  T  (see  Eq.  (1)),  y  will  be  a  strong  function 
of  T.  It  is  therefore  valid  to  define  a  critical  temperature  Tc  that  sepa¬ 
rates  the  two  regimes.  A  rough  estimate  of  Tc  is  obtained  by  setting  the 
non-linear  absorption  coefficient  k,  given  by  Eq.  (1),  to  the  absorption 
coefficient  in  the  inclusion  k^ 

k(Tc)  -  ki 

For  the  purpose  of  definiteness,  in  what  follows  we  shall  set 

AT  -  Tc  -  I*  -  2000  K  . 

If,  as  we  shall  assume,  the  inclusion  has  a  large  thermal  conductivity  as 
compared  to  the  surrounding  medium,  then  the  temperature  profile  during  the 
transient  is  as  shown  in  Figure  3.3.  T  =  Ts  for  r  <  R, 

T  ■  <TS  -  T->4-  *  T-  R  <  r  <  E  . 

and  T  drops  rapidly  to  T«  when  r  >  R  +  /p.  The  thermal  energy  corresponding 
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Figure  3.3.  Temperature  profile  during  heating  stage. 
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to  this  profile  is 


E  -  p^AT  q!*  4nr2dr  +  p^AT  RJ  D  4 nr2  (5] dr 

-  -3-  *Vi"[ 1  *  r  K1  *  V>  2-  ii  ] 


where  the  subscripts  1  and  2  refer  to  the  inclusion  and  to  the  calculus 
material.  Equating  E  with  the  energy  absorbed  from  the  laser  beam,  we  find: 

♦  =  J^Idt'  =  3  RplCl*T  /  l  +  4  [<1  +  *'Xt/R)2  -  l]  ) 

15^1 W~  l  2  piCi  L  J  J 

Ve  set  Qabs  *  1  for  a  strongly  absorbing  inclusion  larger  than  the 
wavelength  and  QabsW  *  (2/3)kjR  for  a  weakly  absorbing  inclusion  (kjR  <  1).‘ 
Ve  show  in  Figures  3.4  and  3.5  the  fluences  for  plasma  ignition  as  a  function 
of  inclusion  radius  for  a  set  of  pulse  lengths.  The  conditions  assumed  for 
the  calculations  were  pj  *  p2  *  2.2  g/cm3,  Cj  *  C2  =  1.5  J  /  (gK),  X  ** 

IQ”2  cm2  s"1  and  AT  -  2000  K. 


1  10  IQ2  103 

LASER  FLUENCE  TO  REACH  aT  »  2000°C  (J/cm2) 


A- 7943 

Figure  3.4.  Laser  fluence  to  reach  AT  =  2000  K, 
strongly  absorbing  particle  (Qabs  =  !)• 
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LASER  FLUENCE  TO  REACH  aT  =  2000° C  (J/cm2) 

A- 7 942 

Figure  3.5.  Laser  fluence  to  reach  AT  =  2000  K, 
weakly  absorbing  particle  (k  *  10^  cm-1). 

One  sees  a  very  different  behavior  of  threshold  with  particle  radius  at 
Tp  *  3  x  10“^s  than  at  the  two  shorter  pulse  lengths  considered.  The  thermal 
diffusion  depth  is,  in  general,  much  larger  (^q  =  20  ym)  than  the  particle 
radius  when  Tp  -  3  x  10-*s,  so  that  the  volume  heated  is  much  larger  than  the 
inclusion.  The  threshold  fluence  will,  under  such  conditions,  vary  inversely 
with  the  power  absorbed  by  the  inclusion.  The  power  absorbed  varies  as 
nR^Qabs,  thus  leading  to  an  R“2  (kR-3)  dependence  for  strongly  (weakly) 
absorbing  inclusions.  If,  however,  <  R,  as  in  the  case  of  the  two  shorter 
pulse  lengths  considered,  the  threshold  fluence  will  vary  as  R  /  Qabs(R)  (see 
Eq.  (10)),  i.e.,  as  R  for  the  strongly  absorbing  inclusions  and  will  be 
independent  of  inclusion  radius  for  the  weakly  absorbing  inclusions.  The 
threshold  irradiance  versus  particle  size  is  shown  in  Figure  3.6  for  a  series 
of  absorption  coefficients. 

It  is  of  interest  to  compare  our  results  with  experimentally  observed 
thresholds  on  calculi.  We  have  done  this  by  introducing  extra  abscissas  in 
Figures  3.4  and  3.5,  corresponding  to  experimentally  observed  breakdown 
fluences  observed  at  the  pulse  lengths  indicated.  All  the  data  were  taken 
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Figure  3.6.  Threshold  irradiance  versus  particle  radius. 


from  the  thesis  of  Watson*.  In  comparing  the  data  vith  the  theory,  one  must 
keep  in  mind  that  experimental  thresholds  probably  include  an  additional 
fluence  for  plasma  growth.  The  threshold  data  on  axis  A  vere  for  uncolored 
calcium  oxalate  calculi  at  A  -  577  nm.  The  threshold  data  on  axis  B,  taken  at 
the  same  wavelength,  were  for  a  dark  pigmented  gallstone.  The  thresholds  for 
B  are  approximately  one  order  of  magnitude  lover  than  for  A. 


Note  that  Teng  et  al  have  also  measured  laser-induced  fragmentation 
thresholds  of  biliary  calculi  at  A  ■  690  nm  and  with  tp  ■  0.8  ns.  From 
Figure  3.1  we  see  that  their  threshold  is  »  4  J/cm^,  vith  plasma  being  formed 
at  t  >  0.4  ps.  If  we  take  into  account  that  their  fiber  diameter  (0.3  mm)  was 
one  half  the  fiber  diameter  of  Watson  and  use  the  threshold  scaling  (derived 
from  the  data  of  Watson)  *  (diameter)-* *3,  ve  obtain  a  threshold  of 
2.4  J/cm^  for  T  «  0.4  ys,  which  is  practically  the  same  threshold  observed  by 
Watson  with  tp  =  1.5  ys.  Because  of  the  slight  wavelength  dependence  of 
fragmentation  threshold  observed  by  Watson,  ve  would  expect  the  threshold  at 
X  -  577  nm  to  be  less  than  at  690  nm  so  that  the  data  of  Teng  et  al  and  Watson 
are  consistent  vith  each  other. 
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Comparison  of  our  theoretical  results  with  experiments  can  be  made  as 
follows.  For  the  strongly  pigmented  gallstone,  the  data  indicate  a  factor  of 
5  increase  in  fluence  threshold  as  one  goes  from  Tp  =  1.5  to  300  us.  Such  a 
spread  in  threshold  is  consistent  with  an  inclusion  radius  of  20  um,  for 
k  >  103  cm-l  (see  Figures  3.5  and  3.4).  The  observed  thresholds,  however,  are 
two  to  four  times  smaller  than  the  calculated  ones.  A  decrease  in  the  temper¬ 
ature  rise  by  a  factor  of  three,  i.e.,  AT  =  700°K  rather  than  AT  =  2000  K,  or 
an  increase  in  k  by  a  factor  of  three  would  bring  the  results  in  better 
agreement.  The  data  for  the  calcium  oxalate  stone  show  practically  the  same 
spread  in  thresholds,  but  threshold  values  that  are  a  factor  of  4  larger  than 
the  model  prediction  for  R  *  20  pm  and  k  =  103  cm-*-.  Agreement  between  the 
model  and  the  data  can  be  achieved  in  this  case  by  decreasing  the  absorption 
coefficient  by  a  factor  of  4.  Though  the  factor  of  ten  increase  in  thresholds 
as  one  goes  from  a  dark  pigmented  gallstone  to  a  clear  calcium  oxalate  urinary 
tract  calculus  can  be  explained  by  a  difference  in  k,  the  expected  lower 
densities  of  absorption  sites  in  the  clear  calculus  should  also  play  a  role. 
Effective  fragmentation  requires  a  pressure  build-up  over  the  whole  irradiated 
surface,  which  requires  that  the  plasmas  ignited  at  the  various  sites  will  be 
larger  in  the  stones  that  have  fewer  ignition  sites. 

We  show  in  Figure  3.7  the  variation  of  breakdown  fluence  with  pulse 
length.  Calculations  were  performed  for  inclusions  of  radius  1,  8  and  20  um. 
The  data  points  have  a  slope  that  is  consistent  with  breakdown  initiated  on 
inclusions  of  radius  20  um,  excepting  for  the  pigment  gallstone  data  point  at 
Tp  »  1.5  us  which  is  a  factor  of  2  too  low. 

We  conclude  that  with  a  proper  (but  reasonable)  choice  of  the  parameters 
in  the  model,  agreement  can  be  achieved  between  the  model  and  the  data.  The 
data  indicate  that  plasma  ignition  is  caused  by  absorbing  inclusions  of  radius 
20  u">  or  larger.  Validation  of  the  model  will  require  a)  better  estimates  of 
the  physical  properties  of  calculi,  b)  evaluation  of  AT,  i.e.,  determination 
of  the  temperature  dependence  of  absorption  coefficient  in  the  bulk  material, 
and  c)  determination  of  the  absorptive  properties  and  size  distribution  of 
inclusions  in  the  various  stones. 


PULSE  TIME  (S) 

A-7944 

Figure  3.7.  Variation  of  breakdown  fluence  with  pulse  length. 
Continuing  Work 

The  rough  calculations  performed  above  indicate  that  plasma  ignition  in 
calculi  may  be  due  to  localized  heating  at  absorbing  sites  of  radius  10-20  urn. 
It  remains  to  be  shown  that  thermal  instabilities  can  lead  to  plasma  growth 
and  overlap  over  the  time  scales  of  the  experiments  in  which  plasmas  and  stone 
shattering  have  been  observed.  Continuing  work  will  involve  solving  the  heat 
conduction  equation  numerically.  A  steady  state  temperature  profile 
corresponding  to  Eq.  (4)  will  be  assumed  initially.  Non-linear  absorption 
with  the  absorption  coefficient  of  the  form  k  -  b  exp  (-E/T)  will  be  assumed 
and  introduced  into  the  heat  source  term.  Estimate  of  the  coefficients  b  and 
E  will  have  to  be  obtained  from  the  literature. 

Analysis  of  Laser  Produced  Plasmas 

Teng  et  al"  have  obtained  spectral  data  of  plasmas  generated  during  laser 
interaction  with  gallstones  and  kidney  stones.  Ve  have  analyzed  these  spectra 
with  the  aim  of  determining  time  evolution  of  plasma  density,  temperature  and 
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pressure.  Ve  reproduce  in  Figure  3.8  Che  spectra  of  plasma  ignited  above 
biliary  calculi  0.22,  0.36  and  1.6  us  after  the  beginning  of  the  laser  pulse. 
The  laser  pulse  duration  was  0.8  us  and  fluence  35  J/cm^.  The  first  two 
spectra  correspond  to  plasma  radiation  during  the  laser  pulse,  while  the  last 
spectrum  vas  taken  0.8  us  after  the  end  of  the  pulse.  Figure  3.8a  shovs  a 
continuum  with  lines  of  Ca  and  Ca+  appearing  in  absorption.  Figures  3.8b  and 
3.8c  also  shov  almost  exclusively  lines  of  Ca  and  Ca+  in  emission.  The 
principal  lines  of  calcium  that  can  radiate  in  the  300-700  nm  range  are  shown 
in  Table  3.1. 

Plasma  density  can  be  obtained  by  looking  at  the  Stark  broadening  of  the 
spectral  lines.  Plasma  temperature  can  be  obtained  by  taking  the  ratio  of 
lines  having  the  same  upper  state  or  by  comparing  the  intensities  of  lines  of 
Ca  and  Ca+  in  order  to  determine  the  ratio  of  ion  to  neutral  populations.  The 
line  strength  for  a  transition  from  an  upper  state  u  to  a  lower  state  1  will 
vary  as  follows: 

Sul*  SUN  Aui  exp  -<EuAT) 

where  Auj  is  the  Einstein  coefficient  for  spontaneous  emission,  gu  the 
degeneracy  of  the  upper  state,  Eu  the  energy  of  the  upper  state,  and  N  the 
number  density  of  radiating  atoms/ions.  The  values  of  g,  A  and  E  for  the 
principal  lines  of  Ca  and  Ca+  are  shown  in  Table  3.1.  The  temperature  can 
also  be  obtained  from  Figure  3.8a  by  comparing  the  spectrum  to  a  black  body, 
assuming  that  the  radiation  is  black  in  the  long  vavelength  part  of  the 
spectrum.  If  ve  can  assume  thermal  equilibrium,  then  the  concentration  of  Ca, 
Ca+,  Ca++  can  be  obtained  by  solving  the  set  of  Saha  equations  for  these 
species.  Note  that  Ca++  is  Argon-like  and  vill  not  produce  line  radiation  in 
the  visible  part  of  the  spectrum.  Ve  shov  in  Figure  3.9  the  degree  of 
ionization  5  of  Ca  as  a  function  of  temperature  for  a  series  of  calcium 
densities.  The  densities  considered  vill  lead  to  electron  concentrations  in 
excess  of  10*®cm“3  and  it  is  important  to  take  into  account  the  lowering  of 
the  ionization  potential  due  to  the  electric  fields  surrounding  the  atoms/ 
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Figure  3.8.  Spectra  of  plasma  radiation  during  and  after  interaction  with  a 
biliary  calculus.  Laser  parameters  tp  *  0.8  ys,  ♦  =  50  J/cm2: 
a)  0.23  ys  after  beginning  of  pulse;  d)  0.36  us  after  beginning 
of  pulse;  c)  1.6  ys  after  beginning  of  pulse  (taken  from 
reference  2). 


TABLE  3.1 

MAIN  LINES  OF  CA  AND  CA+  RADIATING  IN  THE  BAND  300  <  X  <  650  nm 


VAVELENGTH 

(nm) 

Eu 

(cm"1) 

8u 

Aul 

(108  S"1) 

Aul«u 
(10®  S"1) 

Ca 

300.3 

48551 

9 

1.08 

9.7 

363.7 

42746 

15 

0.37 

5.6 

422.6 

23652 

3 

2.18 

6.5 

430.0 

38508 

9 

1.8 

16.0 

444.5 

37754 

15 

0.85 

12.7 

559.2 

38230 

15 

0.56 

8.4 

646.0 

35820 

21 

0.54 

11.34 

Ca+ 

317.3 

56851 

10 

3.62 

36.2 

372.6 

52167 
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2.49 

5.0 
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25414 
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Figure  3.9.  Degree  of  ion  zation  in  calcium  plasma. 


ions.  Ve  used  a  formula  for  the  lowering  based  on  the  Debye  length,  Xq,  in 
the  plasma? 


ap  (j+l)e 

“j 

where  j  is  the  charge  of  the  atom/ion, 


(U) 


Xp  =*  kT  x  ^4n(ne  +  Z  z2  [ca Z+]^' 


-1/2 


and  e  is  the  electron  charge.  One  can  readily  derive  from  z  the  electron 
density  in  a  pure  calcium  plasma  through  the  relation 


n  =  1.5xl022pz  (12) 

In  the  case  of  a  gallstone  or  kidney  stone,  calcium  represents  only  a  fraction 
(f)  per  mass  of  the  material.  Since  Ca  and  Ca*  have  lower  ionization  energies 
(6.1  and  11.9  eV)  than  the  other  atomic  species  (0,  H,  N,  C),  one  expects  that 
calcium  will  be  ionized  well  before  the  others,  so  that,  to  a  first  approxima¬ 
tion,  we  can  still  use  Eq.  (12)  but  must  replace  p  by  fp.  For  a  calcium  oxa¬ 
late  calculus,  the  fraction  per  weight  of  calcium  is  ~20X,  leading  to  f  =0.2. 


Analysis  of  Spectrum  3.8a 

We  estimate  the  electron  density  by  looking  at  the  Stark  broadening  of 
the  two  lines  at  4450A  and  3920A.  The  first  line  is  a  line  of  neutral  calcium 
with  full  width  half  max  of  70A.  From  the  line  broadening  parameters  for  the 
lines  of  Ca  tabulated  by  Griem8,  we  estimate  the  electron  density  to  be 
4  x  1018  cm-3.  At  3920A  we  have  two  lines  of  Ca+  separated  by  35A.  Taking 
into  account  this  separation,  we  estimate  from  -the  spectrum  a  full  width  half 
max  for  each  line  of  65A  corresponding  to  an  electron  density  of  ~10?®  cm~3. 
These  two  densities  are  consistent  with  a  plasma  containing  a  hot  core  that 
radiates  a  continuum,  surrounded  by  a  cooler  layer  containing  mainly  Ca+  and, 
further  out,  mainly  Ca.  The  spectrometer  in  the  experiments  of  Teng  et  al  is 
looking  at  the  plasma  in  a  direction  parallel  to  the  surface.  The  cooler 
layer  corresponds  to  plasma  that  has  expanded  laterally  out  of  the  laser  beam. 
The  electron  density  in  the  heated  inner  core  is  expected  to  be  larger  than 
1020  cm-3.  Ve  have  attempted  to  derive  the  temperature  of  the  inner  core  by 
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assuming  that  radiation  in  the  continuum  is  black.  Ve  have  plotted  in  Figure 
3.10  several  black  body  (BB)  curves  in  units  of  photons  /  (cm^s  nm)  and 
normalized  these  curves  so  that  they  would  pass  through  the  point  having 
ordinate  0.55  at  X  *  450  nm.  Ve  have  also  plotted  in  Figure  3.10  the  spectrum 
corresponding  to  Figure  3.8a.  Ve  see  by  comparing  the  spectrum  to  the  BB 
curves  that  ve  cannot  fit  the  continuum  at  X  >  450  nm  with  any  BB  curve.  The 
continuum  is  inconsistent  with  a  constant  temperature  black  body.  A  possible 
explanation  for  the  more  rapid  decrease  in  photon  flux  with  increasing 
vavelength  observed  as  compared  to  BB  is  the  decreasing  opacity  of  the  plasma 
with  decreasing  vavelength.  The  absorption  coefficient  due  to  free-free  and 
bound-free  IB  collisions  varies  vith  frequency  as^ 

ehvAT 

1 f  S  ,  . . . . 

^7 


The  black  body  flux  in  the  long  wavelength  limit  varies  proportionately  with  T 
for  fixed  v.  The  spectrum  indicates  that  the  effective  temperature  seen  by 
the  radiometer  decreases  by  a  factor  of  at  least  two  as  X  increases  from  450 
to  600  nm. 


Figure  3.10.  Comparison  of  blackbody  spectral  radiation  curves  vith  kidney 
stone  plasma  emission  spectrum  of  Figure  3.8a. 
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The  presence  of  calcium  lines  in  absorption  with  electron  densities  in 
excess  of  10^0  cm-®  is  consistent  with  a  core  density  >  10~2  g/cm®  and  a 
temperature  T  >  15,000  K. 

The  spectrum  in  Figure  3.8b  corresponds  to  a  hotter  but  less  dense  plasma 
than  in  Figure  3.8a.  The  main  lines  of  Ca+  near  3920A  appear  in  emission 
rather  than  in  absorption.  The  feature  at  320  nm  is  an  emission  line  of  Ca+. 
It  is  interesting  to  note  that  this  line  is  absent  from  the  earlier  spectrum 
(Figure  3.8a).  The  spectrum  (Figure  3.8c)  corresponds  to  a  much  cooler  vapor 
since  the  Ca+  line  at  316  nm  has  almost  disappeared.  The  strongest  line  of 
Ca+  at  393  and  396  nm  are  still  apparent,  indicating  that  Ca+  has  not  com¬ 
pletely  recombined.  The  continuum  is  much  weaker  in  the  spectrum,  indicating 
much  lover  electron  densities  than  in  the  first  two  spectra.  This  is  not 
surprising  since  after  the  pulse  the  plasma  has  expanded  and  the  ions  have 
recombined. 

In  order  to  better  analyze  the  spectra,  we  have  generated  spectra  of  pure 
calcium  plasmas  using  the  PSI  radiation  code.  This  code  calculates  emission 
from  isothermal  slabs  of  a  given  thickness.  The  absorption  coefficient  of  the 
plasma  for  free-free,  bound-free,  and  bound-bound  (line)  transition  is 
calculated  and  emission  is  calculated  using  Kirchoff's  lav.  The  formulas  used 
in  the  code  are  described  in  reference  10.  The  data  on  line  strengths  for  Ca 
and  Ca+  was  obtained  from  the  NBS  tabulations1*  and  the  Stark  broadening 
parameters  were  taken  from  the  book  by  Griem.®  We  show  below  in  Figures  3.11 
to  3.16  plasma  radiation  in  the  300-700  nm  band  from  slabs  of  thickness  d. 

The  fiber  used  in  the  experiments  of  Teng  et  al  had  a  diameter  of  0.34  mm.  We 
allowed  the  thickness  d  to  vary  from  0.1  to  1  mm,  which  should  pretty  well 
cover  the  experimental  conditions.  Rather  than  use  thermal  equilibrium 
concentrations,  we  arbitrarily  set  the  concentration  of  Ca,  Ca+  and  Ca++  in 
order  to  better  study  the  effect  of  dominant  ions  on  plasma  radiation. 

Figure  3.11  shows  radiation  from  a  slab  of  thickness  d  >  0.1  mm  at  T  • 
10,000  K  for  a  10%  ionized  calcium  plasma.  We  also  show  in  this  figure  the 
black  body  envelope.  The  slab  is  seen  to  be  optically  thin  except  at  line 
center.  Figure  3.12  shows  the  spectrum  emitted  by  a  plasma  of  the  same 
thickness  as  before  but  with  T  «  25,000  K  and  50%  ionized.  The  concentration 
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Figure  3.11.  Calcium  slab  emission  for  T  »  10,000  K,  p  =  5  Atm, 
ne  -  [Ca+J  «  3.33  x  1017  cm"8  and  [Cal  =  3.3  x  1018 
Slab  thickness  *  0.01  cm. 
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Figure  3.12.  Calcium  slab  emission  for  T  »  25,000  K,  p  =  34  Atm 
ne  -  [Ca+J  -  [CaJ  «  3.33  x  10*8  cm-8.  Slab  thick¬ 
ness  «  0.01  cm. 
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Figure  3.13.  Slab  emission  for  T  *  40,000  K,  Ca+  *  ne  •  1018  cm"3, 
p  =  11  Atm.  Ad  »  0.01  cm,  ad  -  0.03  cm. 
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Figure  3.14.  Slab  emission  for  T  «  40,000  K,  (Ca+J  =  ne 
p  ■  llOAtm.  Ad  =  0.01  cm,  od  =  0.03  cm. 
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Figure  3.15.  Slab  emission  for  T  =  40,000  K,  [Ca+]  =  [ne]  =  1020 
p  *  HOOAtm.  Ad  =  0.01  cm,  ad  =  0.03  cm. 


3000  4000  5000  6000  7000 

WAVELENGTH  (A) 

A-7574 

Figure  3.16.  Slab  emission  for  T  «  40,000  K,  ne  =  [Ca+]  = 
[Ca]  «  3.33  x  1018,  d  -  0.01  cm. 


of  Ca  is  the  same  in  both  spectra.  The  tenfold  increase  in  electron  density 
caused  the  lines  to  broaden  significantly.  Figures  3.13  to  3.16  correspond  to 
a  fully  singly  ionized  calcium  plasma  at  40,000  K.  and  shows  the  effect  of 
density  on  plasma  radiation.  The  Ca+  density  is  allowed  to  vary  from 
10l7  cm-3  (p  =  6.7  x  10-3  g/cm3)  to  10^0  cm-3  (p  =  6.7  x  10“3  g/cm3).  The 
slab  is  seen  to  be  completely  black  at  the  highest  density  considered. 
Introduction  of  an  equal  amount  of  non-ionized  calcium  shown  in  Figure  3.16  is 
seen  to  add  significant  structure  to  the  spectra. 

We  show  in  Figure  3.17  the  photon  flux  emitted  by  a  plasma  slab  of 
density  10~3  g/cm3  and  temperature  15,000  K.  The  composition  chosen  was  that 
corresponding  to  thermodynamic  equilibrium.  The  spectra  of  this  figure  can  be 
directly  compared  with  the  measured  spectra  since  the  detector  response  in 
Teng  et  al's  experiments  was  proportional  to  photon  flux.  We  note  that  the 
Ca+  line  at  320  nm  is  black  and  radiates  more  than  the  fundamental  line  at 
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Figure  3.17.  Photon  flux  emitted  by  a  slab  of  thickness  d  in  thermal 
equilibrium,  T  =  15,000  K,  p  =  60  Atm,  p  =  10-3  g/cm3. 
ne  =  1.41  x  1019,  JT Ca ]  =  1.37  x  1017,  [Ca+]  =  1.33  x  1019, 
[Ca++ ]  =  4.2  x  1017  cm3. 
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400  nm  due  to  the  effect  of  the  black  body  envelope.  The  spectra  in 
Figure  3.8b,  however,  indicates  a  weaker  line  at  320  nm  than  at  400  nm.  The 
only  way  to  reconcile  the  two  spectra  (outside  of  assuming  that  the  spectrum 
of  3.8b  corresponds  to  a  much  lower  temperature,  where  Ca+  should  normally  be 
absent)  is  to  assume  that  the  spectral  response  is  falling  off  at  short 
wavelengths  (due  to  glass  optics).  Teng  et  al,  however,  report  using  quartz 
windows  and  optical  elements. 
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4.  MECHANICAL  FRAGMENTATION  MODELING 


Mechanical  fragmentation  of  kidney  stones  under  laser  irradiation  appears 
to  be  associated  with  the  vaporization  of  a  portion  of  the  calculus  and  the 
subsequent  ionization  and  plasma  formation.  This  high  pressure  plasma  is 
believed  to  generate  a  fragmentation  shock  wave  vhich  propagates  into  the 
stone.  Belov,  a  mathematical  model  is  developed  to  describe  this  phenomena  and 
is  used  to  suggest  the  optimum  laser  waveform  that  will  yield  the  minimum 
energy  requirements  (minimum  energy  per  gram  removed). 

Initially,  laser  pulse  energy  Ep  is  deposited  in  length  over  area  A  to 
obtain  energy  density  e 


A 


where  is  the  effective  optical  absorption  depth,  X,  plus  the  thermal 
diffusion  depth  over  the  pulse  duration  tp, 

=  X  +  vKx 

P 

where  K  is  the  thermal  diffusivity.  When  energy  density  e  exceeds  the  vapori¬ 
zation  threshold,  plasma  formation  is  imminent.  We  write 


po>  is  the  local  bulk  density 
vaporized  by  the  pulse  is 


vhere  Hv  is  the  vaporization  energy  per  gram  and 
of  the  porous  stone  (p®  -  1.5  g/cc).  The  mass  m 
simply  m  =  Ep/Hv. 


From  the  data  of  Watson  (THESIS)  determining  the  fluence  (Fp  =  Ep/A) 
associated  with  the  fragmentation  threshold,  we  deduce 

X  ~  10-20  pm  , 

Hv  -  10  kJ/g  , 
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and 


K  =  0.03  -  0.1  cm2/s 

where  the  value  of  K  is  consistent  with  that  of  0.06  cm^/s  for  limestone  and 
0.02  cm^/s  for  CaO.  The  threshold  fluence  associated  with  vaporization  and 
fragmentation  becomes 


Fth  ■  »JMX  *  *9 

where  we  will  arbitrarily  employ 

Fp  <  l/3Ftjj  :  No  Fragmentation 
l/3Fth  <  Fp  <  Fth  :  Onset  Region 
Fp  *  Fth  :  Fragmentation 

The  actual  pressure  history  and  impulse  delivered  to  the  stone  are 
complicated  functions  of  plasma  formation,  bubble  formation  and  bubble 
collapse.  Vhile  Sections  3  and  5  attempt  to  deal  with  these  issues  in  some 
detail,  a  crude  model  must  be  generated  here  in  order  to  assess  the  shock 
fragmentation  phenomenology. 

If  the  stone  were  irradiated  in  air,  the  vapor  would  leave  the  stone  with 
the  velocity  (V)  of  sound  in  a  gs  vith  the  molecular  weight  of  CaO  at 
approximately  1000°C. 

V  -  4  x  10^  cm/s  . 


The  momentum  Ia  imparted  to  the  stone  is  mV,  and  the  coupling  coefficient  Cj 
becomes 


a 


_  Ia  V  ,  dyne. sec 

CI  *  1!  =  “T  ~  A  Joule 

a  p  v 

and  is  known  as  the  "transparent  vapor  limit"  because  there  is  no  beam 
absorption  above  the  surface  by  either  the  air  or  the  vapor.  Calculations 
appropriate  to  the  limit  of  complete  air/vapor  absorption  (Simons,  1984) 
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suggest  6.9  dyne. sec/ Joule.  In  the  absence  of  any  evidence  of  plasma 
formation  in  air,  the  pressure  on  the  stone  in  air  (pa)  is  determined  in  the 
transparent  vapor  limit 


Noting  that  a  1.5  ys,  30  J/cm^  pulse  satisfies  the  threshold  criterion  (Fth> 
above,  the  pressure  on  the  stone  becomes 


p  =80  bar  . 

*a 

Vhen  the  stone  is  irradiated  under  water,  much  larger  pressures  and 
impulses  are  transferred  to  the  stone.  Teng  et  al  have  observed  that  the 
pressures  under  water  may  be  as  much  as  12  times  greater  than  in  air.  If  the 
corresponding  impulse  coupling  coefficient  is  increased  by  12,  we  obtain 

„  50  dyne . sec 

CI  ~  “Touli -  * 

V 

Until  a  more  rigorous  model  of  bubble  growth  and  collapse  (Section  5)  is 
developed,  the  pressure  imparted  to  a  stone  submerged  in  water  is  assumed  to 
be  constant  during  the  pulse  and  given  by  an  impulse  coupling  coefficient  of 
50  dyne. sec/ Joule. 

»  -  ci  F A 

w 

A  one-dimensional  slug  model  is  used  to  describe  the  interaction  of  a 
high-intensity  laser  pulse  with  a  porous  kidney  stone.  The  mass  loss  due  to 
vaporization  is  neglected  and  the  interaction  is  approximated  by  replacing  the 
laser  pulse  by  a  one-dimensional  piston  pushing  on  the  target  surface  with  the 
equivalent  pressure.  Vhen  the  pressure  is  greater  than  the  strength  of  the 
porous  stone,  ou,  a  fracturing  shock  wave  will  propagate  into  the  target.  The 
zone  between  the  shock  front  and  the  piston  is  filled  with  the  fragmented 
target  debris.  The  analytic  model  (Simons  and  Legner,  1982)  for  the  Shock 
Hugoniot  in  a  porous  material  is  used  to  obtain  the  conditions  behind  the 
shock.  The  density  p  of  the  partially  compressed  fragmented  porous  material 
behind  the  shock  is  given  by 
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p  l1  +  j]  pKo  +  1 

where  p  is  the  pressure  induced  by  the  laser  pulse ,  K0  is  the  cold  compression 
coefficient  of  the  kidney  stone,  r  is  the  GrQneisen  coefficent,  6  is  the 
porosity,  and  p0  is  the  non-porous  density  of  the  stone.  The  fracture  shock 
velocity,  Vf,  can  then  be  written  as 


V 


f 


r  PP  11/2 

Lp>  -  pJJ 


where  pm  is  the  bulk  density  of  the  porous  material  in  front  of  the  shock. 

PaD  -  po(l-0)  . 

Since  the  physical  properties  of  kidney  stones  are  not  readily  available, 
the  properties  of  CaO  will  be  used  whenever  possible.  Hence,  we  use 


pQ  -  3g/cc 
0  -  0.5 
pm  -  1.5g/cc 

r  -  2* 

K0  ~  3.5  x  10“12  cm^/dyne 

and, 


<xu  «  80  bar  . 


where  the  choice  of  80  bar  for  ou  is  non-trivial.  Handbook  values  for  the 
compressive  strength  of  concrete  and  limestone  are  upwards  of  350  bar  with 
corresonding  values  of  the  limiting  shear  at  80  bar.  Since  kidney  stones  are 
approximately  50%  porous,  we  consider  their  ultimate  strengths  to  be  175  bar 
in  compression  and  40  bar  in  shear.  Vhen  a  material  is  subjected  to  a 
one-dimensional  compression  of  strength  a,  the  corresponding  shear  45°  off 
axis  is  c/2.  Hence,  when  o  is  80  bar,  the  material  will  fail  in  shear  at  45° 


*  varies  between  1.5  and  2.5  for  metals  and  reduces  to  r-1  for  a  gas.  Ve 
assume  a  value  of  2  for  f. 
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to  the  principal  load  and  not  in  compression  on  the  principal  axis.  This 
value  of  ou  is  supported  by  compressive  tests  we  have  performed  on  stones  in 
our  laboratory,  indicating  failure  in  the  range  of  60  bar  to  110  bar. 

At  the  end  of  the  laser  pulse,  a  high  pressure  slug  of  fragmented  target 
material  of  density  p  is  lodged  in  the  stone  as  illustrated  in  Figure  4.1. 
Following  the  deposition  of  energy  by  the  laser,  the  pressure  constraint 
(i.e.,  the  piston)  at  the  interface  is  removed  and  a  rarefaction  wave  begins 
to  propagate  toward  the  shock  front.  The  shock  wave  continues  to  propagate 
into  the  target  until  the  rarefaction  wave  reaches  the  shock  front  and  reduces 
its  pressure  to  ou. 
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Figure  4.1.  Stone  fragmentation  at  the  end  of  the  laser  pulse. 
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For  the  description  of  the  flovfield  after  the  laser  pulse  is  off,  a 
simplified  Arbitrary  Lagrangian-Eulerian  computer  program  (Amsden  et  al., 
1980),  SALE,  is  used.  The  hydrodynamic  code-SALE  was  originally  developed  at 
the  Los  Alamos  Scientific  Laboratory  and  is  available  at  the  National  Energy 
Software  Center,  Argonne  National  Laboratory,  Argonne,  Illinois.  The  SALE 
code  combines  an  implicit  treatment  of  the  pressure  equation  similar  to  that 
in  the  Implicit  Continuous-Fluid  Eulerian  (ICE)  technique  with  the  grid 
rezoning  philosophy  of  the  Arbitrary  Lagrangian-Eulerian  method.  As  a  result, 
it  can  handle  flow  speeds  from  supersonic  to  the  incompressible  subsonic  limit 
in  a  grid  that  may  be  moved  with  the  fluid  in  typical  Lagrangian  fashion,  held 
fixed  in  a  Eulerian  manner,  or  moved  in  some  arbitrary  way  to  give  a 
continuous  rezoning  capability. 

The  SALE  code  offers  flexibility  in  specifying  the  equation  of  state  of 
the  shock  layer.  For  the  present  calculations,  an  ideal  gas  assumption  is 
made  for  the  fragmented  target  material.  In  general,  the  "SESAME  TABLE,” 
which  is  the  tabulated  equation  of  state  for  numerous  materials  compiled  by 
the  Los  Alamos  National  Laboratory,  may  be  used. 

Illustrative  examples  of  the  results  of  the  SALE  code  have  been  obtained 
for  the  sample  case  of  a  1  iis  pulse  with  a  fluence  of  200  J/cm^,  irradiated 
onto  a  50  percent  porous  kidney  stone.  The  density,  pressure,  and  velocity 
contours  corresonding  to  a  time  0.64  ys  after  the  pulse  are  illustrated  in 
Figures  4.2  through  4.4,  respectively.  The  contours  clearly  demonstrate  that 
the  fragmented  target  debris  has  expanded  well  outside  the  crater.  This 
expansion  enhances  the  dissipation  of  the  impact  energy  through  the  debris 
velocity. 

The  goal  of  the  present  calculations  is  to  assess  the  effects  of  the 
pulse  fluence  and  pulse  time  on  the  fracturing  of  kidney  stones  by  laser 
radiation.  Figure  4.5  illustrates  the  crater  depth  as  a  function  of  pulse 
time  for  various  pulse  fluences.  For  a  pulse  fluence  of  200  J/cm^,  the  crater 
depth  increases  from  about  4  x  10“2  cm  at  Tp  =  10"^s  to  approximately  1.5  cm 
at  Tp  *  10~*s.  Since  the  pulse  fluence  and  the  coupling  coefficient  are 
assumed  to  be  constant,  the  pressure  decreases  with  increasing  pulse  duration 
(decreasing  laser  intensity).  In  the  one-dimensional  limit,  there  is  no 
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Figure  4.2.  Density  contours  at  1  us  after  laser  energy 
deposition. 
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Figure  4.4.  Axial  velocity  contours  at  1  ys  after  laser  energy 
deposition. 
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Figure  4.5.  One-dimensional  crater  depth. 
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geometric  decay  of  the  fragmentation  shock  and  lower  pressures  applied  over 
longer  times  yield  deeper  craters.  The  deepest  crater  then  occurs  when  the 
induced  pressure  just  exceeds  the  strength  of  the  stone.  However,  the  three- 
dimensional  decay  of  the  pressure  pulse  will  seriously  alter  this  result.  We 
presently  plan  to  treat  three-dimensional  effects  in  the  coming  year. 

The  present  calculations  also  demonstrate  the  obvious  result  that  the 
crater  depth  decreases  with  decreasing  pulse  fluence  until  that  fluence  is 
below  the  fragmentation  threshold.  The  dashed  lines  represent  calculations  in 
the  "onset  region" 

Fth  <  Fp  <  Fth  . 

The  energy  required  for  the  mass  removal  is  defined  as 

F 

E  =  P 
t)  p  x  Depth 

and  is  shown  in  Figure  4.6.  For  a  pulse  fluence  of  200  J/cm2,  the  energy 
required  to  remove  a  unit  mass  of  stone  decreases  from  3  x  10^  J/gm  at  10_7s 
to  about  102  J/gm  at  10~4s.  The  minimum  energy  requirement  occurs  at  the 
point  where  the  pressure  provided  by  the  laser  pulse  is  just  above  threshold. 
This  is  largely  due  to  the  one-dimensional  approximation.  The  energy 
efficiency  also  increases  (energy  decreases)  with  decreasing  pulse  fluence 
until  below  threshold.  Again,  the  dashed  lines  represent  the  pulse  fluence  in 
the  range 

-3  F<h  <  F<>  <  F,h  ' 

The  above  calculations  suggest  a  maximum  energy  efficiency  (minimum 
energy/mass)  when  operating  just  above  the  threshold  of  mechanical  fragmen¬ 
tation.  These  predictions  are,  however,  subject  to  many  approximations: 

1)  one-dimensional  geometry 

2)  constant  pressure  over  the  pulse  duration 

3)  50  dyne. sec/ Joule  coupling  coefficient 

4)  porous  Shock  Hugoniot  for  porous  volumes  that  are  filled  with  air. 

5)  the  neglect  of  spall  phenomena 
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ENERGY  REQUIRED  FOR  MASS  REMOVAL  (J/GM) 


5.  A  SIMPLIFIED  MODEL  FOR  LASER  LITHOTRIPSY 
-  Modeling  of  the  Above-Surface  Fluid  Dynamics  - 

Laser  lithotripsy  appears  to  depend  upon  laser-induced  breakdovn  at  or 
near  the  surface  of  a  high  pressure  bubble.  The  resulting  breakup  of  the 
stone  depends  upon  the  resulting  spatial  and  temporal  pressure  history  which 
the  stone  is  subjected  to.  The  purpose  of  this  analysis  is  to  suggest  a 
simplified  method  for  treating  the  fluid  mechanics  of  the  problem.  Hopefully, 
this  would  be  helpful  both  in  optimizing  laser  pulse  characteristics  and  in 
using  acoustic  signals  to  assist  in  assessing  effectiveness. 

The  overall  flow  problem  may  be  separated  into  treatment  of  the  vapor 
bubble  on  the  one  hand  and  treatment  of  the  resulting  flow  in  the  essentially 
incompressible  fluid  surrounding  it  on  the  other.  Since,  in  general,  the 
fluid,  which  we  will  assume  to  have  the  properties  of  water,  represents 
significant  inertia,  its  flow  over  a  wide  range  of  conditions  can  be 
approximated  quite  well  by  a  linearized  or  acoustic  model.  Thus,  we  will 
assume  that  the  water  flow  is  described  by  sound  waves  spreading  outward  from 
the  vapor  bubble. 

The  boundary  condition  at  the  surface  of  the  bubble  is  that  the 
pressure  and  expansion  velocity  of  the  bubble  must  match  the  pressure  and 
fluid  velocity  of  the  water. 

The  vapor  bubble,  on  the  other  hand,  will  tend  to  expand  slowly 
compared  to  the  sound  speed  in  the  vapor.  Thus,  the  vapor  bubble  may  be 
treated  as  a  uniform  pressure  region.  The  properties  of  the  vapor  will  be 
determined  by  absorption  of  laser  energy  and  by  loss  of  energy  by  reradiation 
and  by  PV  work  done  on  the  surrounding  water. 

The  assumptions  introduced  thus  far  are  valid  if  the  expansion 
velocity  of  the  bubble  is  small  compared  to  the  sound  speed  in  both  the  water 
and  the  vapor  bubble.  Actually,  the  results  may  be  somewhat  forgiving  as  to 
violation  of  the  first  condition  since  the  radial  expansion  of  the  wave  in  the 
vater  will  rapidly  reduce  the  strength  of  the  wave. 
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Vapor  Bubble 


For  the  present  analysis,  ve  will  further  assume  that  the  temperature 
and  density  within  the  vapor  bubble  are  uniform  and  that  the  heat  capacity  can 
be  represented  by  a  constant.  The  latter  assumption  neglects  heats  of 
ionization,  dissociation,  and  vaporization. 


The  energy  equation  is  then  given  by 

d< PcvTV)  dV  dE 

- 3t -  +  p  3t  "  Si 


(1) 


or 


r 

r^T 


+ 


1 


dE 

*  Jt 


(2) 


where  p  is  the  density,  Cy  is  the  heat  capacity  at  constant  volume,  T  is  the 
temperature,  V  is  the  total  volume  of  the  gas  bubble,  p  is  the  pressure,  dE/dt 
is  the  rate  of  absorption  of  laser  energy  (or  minus  the  rate  of  reradiation), 
and  y  is  the  ratio  of  specific  heats. 


For  the  case  of  dE/dt  »  0,  this  reduces  as  expected  to 


pVY  -  constant  *  p  V  Y  (3) 

o  o 

Note  that  within  the  approximation  of  negligible  heat  of  vaporization  this  is 
true  independent  of  the  amount  of  vaporization  which  has  occurred. 

One-Dimensional  Planar  Geometry 

For  illustrative  purposes,  let  us  first  consider  the  planar  one¬ 
dimensional  case.  In  this  case,  V  and  E  are  to  be  interpreted  as  the  volume 
and  net  energy  absorbed  per  unit  area  respectively.  Thus,  V  becomes  simply 
the  linear  thickness  of  the  vapor  bubble  and  the  velocity  at  the  edge  of  the 
bubble  is  given  by 
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u  =  dV/dt  (4) 

For  waves  moving  outvard  only  (towards  positive  x),  the  pressure  and 
velocity  in  the  water  are  functions  of  [ct  -  x]  only  and  not  of  [ct  +  x], 
where  c,  t,  and  x  are  the  speed  of  sound  in  the  water,  the  time,  and  the 
position  respectively.  Further,  the  acoustic  equations  give 

5p  =  p^cSu  (5) 

as  the  change  across  each  elemental  wave  (where  is  now  the  density  of  the 
liquid).  Taking  the  initial  uniform  pressure  as  pa  and  the  initial  velocity 
as  zero  throughout  the  medium,  we  may  integrate  (5)  to 

p  =  ^cu  +  p^  (6) 

Vi  thin  the  acoustic  approximation,  the  vapor-fluid  boundary  is 
essentially  at  x  =  0.  Using  equations  2,  4,  and  6  and  matching  velocity  and 
pressure  at  the  boundary,  we  obtain 


dE 

Jt 


(7) 


for  the  vapor  volume.  Since  the  above  relations  relate  V  and  its  derivatives 
to  pressure  and  velocity,  and  since  these  are  constant  along  lines  of  constant 
ct-x,  integrating  Eq.  (7)  for  a  given  dE/dt  defines  the  entire  flow  field. 
Thus,  provided  that  a  reasonable  approximation  to  dE/dt  exists,  numerical 
integration  of  Eq.  (7)  will  provide  a  solution. 


For  the  present,  let  us  consider  only  the  highly  simplified  case  of  a 
very  short  laser  pulse  and  examine  the  flow  field  only  in  the  period  after  the 
end  of  light  absorption,  dE/dt  =  0.  In  this  case,  we  may  use  equations  3,  4, 
and  6  to  obtain 

(■>-  *  ri?)vT  ■  "oV  <8> 

or 

f  ar  =  po(~r)  "  p®  (9) 
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Thus,  when  the  pressure  is  initially  large  compared  to  p<»,  the  volume  will 
grow  as  tlA'+l  and  the  pressure  will  decay  as  t~YA’+l.  At  late  times,  the 
volume  will  asymptotically  approach  ^0(p0/pa)^^r,  corresponding  to  the 
pressure  asymptotically  approaching  the  ambient  pressure  pa. 

Experimental  observations  of  roughly  spherical  bubble  formation  indicate 
that  the  bubble  reaches  a  maximum  size  and  then  collapses.  In  the  present 
planar  case,  the  theory  does  not  predict  the  collapse.  The  difference  is 
presumably  due  to  the  difference  in  the  number  of  dimensions. 


Spherically  Symmetric 


In  the  spherically  symmetric  case,  the  acoustic  solution  for  outgoing 
waves  only  is  given  by 

r+  -  f[ct-r]  (10) 

where  f  is  an  arbitrary  function  of  (ct-rj,  r  is  the  radial  coordinate, 

94  f ' 

P  *  Pfgj  +  const  »  ye-  +  pm  (11) 


and 


(12) 


At  the  edge  of  the  bubble,  r  =  a,  we  may  assume  that  f'/a  «  f/a^  since  f  will 
vary  on  a  scale  of  ct,  which  under  the  acoustic  assumption  is  large  compared 
to  a.  Thus,  the  velocity  at  the  edge  of  the  bubble,  Ue,  can  be  approximated  as 


Since 


(13) 


(14) 


(15) 


Eq.  (13)  becomes 

,  1  dV 

£  =  ~TR  ST 

and  using  Eq.  (11),  Eq.  (2)  becomes 


1 

FT 


(r-1/3  )Po„- 1/3 

- Z'T/3*  v 

(48jt  ;  J 


d2v 


dt 


+  YP0 


dV 

ar 


.2/3  d3V 


IF 


dE 
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(16) 


which  again  can  be  integrated  numerically  for  a  given  dE/dt  to  define  the 
entire  flow  and  pressure  field. 


If  we  look  again  at  the  simplified  case  of  a  very  short  laser  pulse  so 
that  in  the  solution  region  dE/dt  -  0,  Eq.  (3)  becomes 


d2V 


- 7?  2/3 

(48lt  )  dt 


(17) 


Multiplying  both  sides  by  dV/dt  and  integrating,  one  obtains 


+  A 


(18) 


The  integration  constant,  A,  would  be  determined  by  the  initial  conditions, 
i.e.,  integration  of  Eq.  (16)  during  the  laser  pulse.  Thus  far,  we  have  been 
unable  to  find  a  satisfactory  approximation,  short  of  this  integration,  for 
determining  A. 


Prom  the  form  of  Eq.  (18)  it  is  possible  to  conclude  that,  since  for 
realistic  values  of  y,  y-4/3  is  a  small  number,  dV/dt  will  initially  be  a 
constant.  Also,  since  it  is  the  square  of  dV/dt  that  approaches  zero  at  some 
point,  the  solution  will  show  a  maximum  volume  followed  by  a  collapse.  Note 
that  the  maximum  volume  is  not  determined  by  a  simple  relation  of  pressure  to 
p„  but  is  dependent  on  initial  conditions. 
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Equation  18  seems  to  shov  a  rather  strong  dependence  on  initial 
conditions  and  is  not  determined  solely  by  the  total  laser  energy  absorbed. 
It  is  tempting  to  speculate  that  the  pressure,  as  veil  as  its  spatial  and 
temporal  variation,  and  therefore  the  effectiveness  of  the  laser  lithotripsy 
may  be  strongly  dependent  on  detailed  characteristics  of  the  laser  pulse. 
Certainly  it  would  seem  that  integration  of  Eq.  (16)  for  a  variety  of  energy 
absorption  curves  would  be  useful. 


Discussion 


The  above  analysis  has  not  discussed  the  mechanisms  of  laser  light 
absorption.  Clearly,  a  more  sophisticated  consideration  of  vapor  properties 
similar  to  that  used  in  other  laser  materials  interaction  studies  is  required. 
The  major  difference  in  consideration  of  vapor  properties  is  that  in  this  case 
we  are  suggesting  that  the  vapor  expands  slowly  under  quasi-hydrostatic 
conditions  rather  than  the  sonic  expansion  velocities  which  are  reached  when 
the  laser  material  interaction  occurs  without  a  surrounding  high  density 
material. 


The  use  of  the  acoustic  approximation  for  the  flow  in  the  surrounding 
medium  may  be  helpful  in  defining  the  usefulness  of  pressure  monitors  placed 
at  some  distance  from  the  sight  of  the  laser  focus  for  measuring  the  strength 
of  the  laser-induced  pressure  pulse.  This  could,  of  course,  be  limited  by 
non-uniform  acoustic  propagation  and  reflections. 

We  do,  however,  hope  that  the  basic  approximation  of  separating  the 
flow  into  a  quasi-static  vapor  bubble  and  an  acoustic  flow  field  will  prove 
useful  as  each  of  the  separate  regions  is  considered  in  more  detail. 
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